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Preface 


Current  numeric  calculations  of  state-to-state  resolved  quantum  reactive  scat¬ 
tering  matrix  elements  are  computationally  limited  to  systems  of  three  degrees  of 
freedom  or  less.  The  large  grids  or  a  large  number  of  grids  used  in  the  computation 
are  inefficient.  In  an  effort  to  improve  the  computation  times  in  two-  dimensional 
calculations,  the  combination  of  the  channel  packet  method  together  with  absorb¬ 
ing  boundary  conditions  is  explored  for  the  collinear  H  +  H2  reaction  and  a  simple 
model  of  two  collinear  coupled  Morse  oscillators.  For  collinear  H  +  H2,  results  are  ob¬ 
tained  which  are  in  good  agreement  with  previous  calculations  while  simultaneously 
providing  a  reduction  in  computation  time.  For  the  simple  model  of  two  coupled 
collinear  Morse  oscillators,  scattering  matrix  elements  for  some  mass  configurations 
contain  non-trivial  errors  introduced  by  absorbing  boundary  condition  reflection. 
For  these  mass  configurations,  the  effects  of  absorbing  boundary  condition  reflection 
are  explored  and  scattering  matrix  elements  are  presented  for  absorbing  boundary 
condition  reflection  for  a  simplified  Morse  potential.  Good  scattering  matrix  ele¬ 
ments  are  presented  for  the  light-heavy-light  mass  configuration  of  the  two  coupled 
collinear  Morse  oscillators.  The  effects  of  kinetic  coupling  and  potential  well  depth 
on  scattering  matrix  elements  are  also  presented  for  this  mass  configuration. 
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Abstract 

In  an  effort  to  develop  a  more  efficient  time  dependent  approach  for  calculating 
scattering  matrix  elements,  absorbing  boundary  conditions  are  combined  together 
with  the  channel  packet  method.  The  channel  packet  method  relies  on  calculat¬ 
ing  two  M0ller  states:  one  representing  reactants  and  one  representing  products. 
The  time  dependent  correlation  function  between  the  two  Mpller  states  is  efficiently 
computed  by  individually  propagating  the  Mpller  states  using  absorbing  boundary 
conditions  as  they  exit  the  interaction  region  of  the  potential.  As  a  Mpller  state 
evolves  in  time,  it  will  be  attenuated  by  the  absorbing  boundary  conditions  without 
reflecting  off  the  edge  of  the  grid  thereby  permitting  the  use  of  a  much  smaller  grid. 
The  Fourier  transform  of  the  correlation  function  is  then  used  to  compute  scattering 
matrix  elements. 

The  efficiency  of  the  combination  of  the  channel  packet  method  with  absorbing 
boundary  conditions  is  demonstrated  in  one  dimension  through  application  to  the  one 
dimensional  square  well.  The  one  dimensional  square  well  has  an  analytic  solution 
which  provides  a  benchmark  for  testing  both  the  efficiency  and  accuracy  of  this 
approach.  Scattering  matrix  elements  obtained  using  the  channel  packet  method 
with  absorbing  boundary  conditions  are  in  excellent  agreement  with  the  analytic 
solution.  The  new  approach  converges  to  the  correct,  analytic  solution  cis  well  as 
provides  a  dramatic  reduction  in  grid  size.  The  reduced  grid  size  results  in  a  faster, 
more  efficient  computation.  The  effects  of  absorbing  boundary  conditions  in  one 
dimension  are  investigated  using  a  one  dimensional  potential  consisting  of  a  Gaussian 
well  with  symmetric  Gaussian  barriers. 

Similar  to  the  one  dimensional  square  well,  the  collinear  H  H2  reaction 
provides  a  benchmark  for  testing  the  efficiency  and  convergence  of  the  channel  packet 
method  with  absorbing  boundary  conditions  as  applied  to  scattering  matrix  elements 

xxii 


in  two  dimensional  scattering.  Again,  the  new  method  is  more  efficient  than  previous 
methods  while  converging  on  the  accepted  solution.  The  combination  of  the  channel 
packet  method  with  absorbing  boundary  conditions  results  in  an  order  of  magnitude 
savings  in  the  time  necessary  to  compute  the  correlation  function. 

With  the  channel  packet  method  with  absorbing  boundary  conditions  estab¬ 
lished  as  an  efficient,  accurate  method  for  computing  quantum  reactive  scattering 
matrix  elements,  a  model  system  consisting  of  two  coupled  collinear  Morse  oscilla¬ 
tors  is  used  to  investigate  the  effects  of  kinetic  energy  coupling  and  potential  well 
depth  on  scattering.  However,  for  light-light-light,  medium-light-medium  and  heavy- 
light-heavy  mass  configurations,  absorbing  boundary  condition  reflection  introduces 
significant  error  into  the  scattering  matrix  elements.  In  an  effort  to  understand 
the  errors  introduced  by  absorbing  boundary  condition  reflection,  scattering  matrix 
elements  are  presented  for  reflection  from  absorbing  boundary  conditions  using  a 
simplified  Morse  oscillator  potential.  Unlike  the  previous  three  mass  configurations, 
scattering  matrix  elements  for  the  light-heavy-light  mass  configuration  do  not  suf¬ 
fer  from  significant  absorbing  boundary  condition  reflection  error.  For  this  mass 
configuration,  the  effects  of  kinetic  energy  coupling  and  well  depth  on  scattering 
matrix  elements  are  presented  for  a  variety  of  two  collinear  coupled  Morse  oscillator 
potentials. 
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REACTIVE  QUANTUM  SCATTERING  IN  TWO  DIMENSIONS 


/.  Introduction 

Molecular  reaction  dynamics  plays  a  central  role  in  many  research  efforts 
of  interest  to  the  United  States  Air  Force.  For  example,  the  High  Energy  Density 
Materials^’^  project  at  Edwards  AFB,  CA  seeks  to  identify  new  rocket  fuels  that  will 
improve  the  current  spacelift  capabilities  of  the  Air  Force.  This  effort  is  supported 
and  guided  by  a  variety  of  computational  models  that  primarily  focus  on  determining 
the  molecular  structure  and  energies  of  high  energy  density  material  candidates. 
Both  classical  and  quantum  mechanical  models  of  the  dynamical  properties  will 
play  an  essential  role  in  determining  the  suitability  of  a  new  high  energy  density 
material  for  use  as  a  rocket  fuel.  Another  research  effort  of  interest  to  the  Air  Force 
is  the  development  of  more  powerful  and  compact  chemical  lasers  at  the  Air  Force 
Institute  of  Technology  and  Phillips  Lab,  Kirtland  AFB,  Chemical  lasers 

will  be  useful  in  a  variety  of  Air  Force  applications  including  the  Airborne  Laser 
project  and  the  Infrared  Countermeasures  project.^”^  The  continued  development  of 
chemical  lasers  will  benefit  from  a  more  detailed  view  of  the  collisional  transfer  of 
energy  afforded  by  molecular  dynamics  simulations.  In  addition,  several  Air  Force 
research  projects  at  the  Geophysics  Directorate,  Phillips  Lab,  Hanscom  AFB,  MA 
are  investigating  the  dynamics  of  upper  atmosphere  pollution.®”*®  The  research  into 
upper  atmosphere  pollution  is  concerned  with  both  the  sources  of  pollution  as  well 
as  the  many  chemical  reactions  between  pollutants  and  the  natural  constituents  of 
the  upper  atmosphere.  For  example,  what  is  the  impact  of  high  altitude  flights  by 
supersonic  aircraft,  military  or  commercial,  on  the  ozone  layer?  The  answer  to  this 
question  and  many  others  will  benefit  from  the  ability  to  efficiently  model  reaction 
dynamics  involved  in  upper  atmosphere  pollution. 
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The  quantum  dynamics  of  a  molecular  reaction  can  be  understood  through 
the  calculation  of  quantum  scattering  matrix  (S-matrix)  elements.  The  absolute 
value  squared  of  an  S-matrix  element  is  the  state-to-state  resolved  probability  of 
reaction.  Current  methods  for  computing  S-matrix  elements  are  limited  to  very 
simple  reactions  where  the  number  of  dimensions  is  two  or  three.  The  research 
presented  here  is  part  of  an  effort  to  develop  computational  methods  which  are  more 
efficient  than  current  methods  and  thereby  extend  the  range  of  reactions  for  which 
S-matrix  elements  may  be  calculated. 

1.1  Background 

It  is  only  in  the  last  two  decades  that  time  dependent  methods  have  become 
practical  for  computing  S-matrix  elements.  Historically,  the  first  practical  methods 
for  computing  S-matrix  elements  were  analytical  solutions  to  very  simple  problems. 
In  1929,  London  first  presented  his  method  for  analytically  modeling  potential  energy 
surfaces. London’s  method  only  took  into  account  three,  collinear  theoretical  atoms 
A,  B  and  C  and  yields  an  analytic  solution  for  the  S-matrix  elements.  In  1931, 
Eyring  and  Polanyi  e.xtended  London’s  methods  to  produce  an  analytical  surface 
and  S-matrix  elements  for  the  two  dimensional,  collinear  H  +  H2  reaction. The 
LEP  surface,  as  it  was  named,  correctly  matched  several  features  of  the  true  surface; 
though,  the  surface  incorrectly  predicts  a  potential  well  for  the  intersection  of  H2  +  H 
and  H  A  H2.  It  wasn’t  until  1955  when  Sato^^  corrected  the  LEP  surface  (to  the 
LEPS  surface)  to  show  a  barrier  in  the  interaction  region  of  H  AH2  —^H2AH.  The 
analytical  methods  of  London,  Eyring,  Polanyi  and  Sato  remained  state  of  the  art 
for  computing  S-matrix  elements  well  into  the  1970s. 

In  the  1950’s,  new  methods  that  relied  on  the  numeric  solution  of  time  in¬ 
dependent  classical  and  quasi-classical  equations  became  possible  with  the  advent 
of  the  digital  computer. The  disadvantage  to  these  methods  is  that  they  rely 
upon  the  diagonalization  of  large  matrices  whose  elements  are  the  coefficients  of  the 
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equations  of  motion  of  the  reacting  particles.  The  disadvantages  of  large  grids  and 
matrices  with  the  high  financial  cost  of  computer  time  continued  to  limit  their  use  to 
simple  one  and  two  dimensional  problems.  In  the  early  1980s  faster,  less  expensive 
computers  became  available  and  quantum  S-matrix  calculations  became  computa¬ 
tionally  feasible. However,  these  calculations  are  limited  in  scope  to  either  two 
dimensional  geometries  for  a  variety  of  reactions  or  three-dimensional  calculations 
for  very  simple  reactions  like  D  +  H2  DH  -f-  H.  While  some  quantum  meth¬ 
ods  allow  the  computation  of  a  single  state- to-state“  resolved  S-matrix  element,  the 
computational  overhead  associated  with  the  large  grids  necessary  to  perform  these 
calculations  remains  the  limiting  factor.  A  reduction  in  grid  size  while  maintaining 
the  accuracy  of  the  result  will  lead  to  faster  computations  and  the  ability  to  compute 
S-matrix  elements  for  a  wider  variety  of  reactions. 

1.2  Problem  and  scope 

The  goal  of  this  research  is  to  develop  the  combination  of  the  time  depen¬ 
dent  channel  packet  method  (usually  referred  to  as  the  channel  packet  method)^^’^'’^® 
with  absorbing  boundary  conditions'’®’ as  a  valid,  accurate  and  efficient  approach 
to  calculating  quantum  scattering  matrix  elements  in  one  and  two  dimensions.  The 
new  method  will  be  tested  for  validity  through  the  application  of  several  conver¬ 
gence  tests.  The  first  tests  convergence  towards  the  correct  solution.  A  method 
which  converges  to  a  wrong  answer  has  no  use.  In  one  dimension,  the  square  well 
potential  has  an  analytic  solution  a  provides  an  excellent  test  for  convergence  to  the 
correct  solution.  In  two  dimensions,  the  collinear  H  -{■  H2  reaction  is  well  estab- 
lished’'’^®’^’”®®’^®  as  a  benchmark  for  computational  quantum  chemistry  and  will  be 
used  to  test  the  new  method  for  convergence  towards  the  correct  solution.  The  test 
for  convergence  to  the  correct  solution  also  demonstrates  the  accuracy  of  the  new 
method.  The  second  test  for  convergence  is  concerned  with  efficiency,  or  order.  A 

“Where  the  states  involved  are  limited  to  those  described  in  the  Hamiltonian. 
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method  can  be  valid  while  also  being  cumbersome  and  computationally  expensive. 
The  application  of  the  channel  packet  method  with  absorbing  boundary  conditions 
will  be  tested  for  efficient  convergence  for  the  square  well  potential  in  one  dimension 
and  the  collinear  H  +  H2  potential  in  two  dimensions.  The  efficiency  of  the  new 
method,  gained  by  reducing  the  size  of  the  grids  used  in  the  computation,  can  be 
judged  by  the  time  a  computation  using  the  channel  packet  method  takes  with  and 
without  absorbing  boundary  conditions.  Finally,  the  influence  of  absorbing  bound¬ 
ary  conditions  on  the  computational  results  will  be  investigated  using  a  Gaussian 
potential  in  one  dimension  and  two  coupled  Morse  oscillators  in  two  dimensions. 

The  computer  resources  used  for  this  study  include  the  MIPS-R3000-based 
Silicon  Graphics  workstations  at  the  Air  Force  Institute  of  Technology,  the  MIPS- 
RlOOOO-based  Silicon  Graphics  Power  Challenge  Array  at  the  Major  Shared  Resource 
Center  at  Wright-Patterson  AFB  and  the  MIPS-R5000-based  Silicon  Graphics  work¬ 
stations  at  the  Air  Force  Research  Lab,  Hanscom  AFB.  The  Fortr.\n77  computer 
language  was  used  for  all  code  development. 

1.3  Outline 

Chapter  II  provides  a  basic  overview  of  quantum  reactive  scattering,  ab¬ 
sorbing  boundary  conditions,  the  channel  packet  method  and  the  underlying  use  of 
Mpller  operators,  and  the  numeric  implementation  of  the  channel  packet  method 
with  absorbing  boundary  conditions.  Chapter  III  presents  the  results  of  scattering 
in  one  dimension  for  a  square  well  potential  and  a  Gaussian  barrier-well-barrier  po¬ 
tential.  Chapter  IV  researches  the  application  of  the  combination  of  the  channel 
packet  method  with  absorbing  boundary  conditions  to  two  dimensional  quantum 
scattering  by  computing  S-matrix  elements  for  the  collinear  H  +  H2  reaction.^* 
The  problem  of  reflection  from  problems  of  absorbing  boundary  conditions  is  re¬ 
searched  in  Chapter  V  for  a  model  two  dimensional  system  of  two  coupled  Morse 
oscillators.  Absorbing  boundary  condition  reflection  introduces  error  into  the  S- 
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matrix  elements  for  certain  mass  configurations.  As  an  initial  analysis  of  reflection 
from  absorbing  boundary  conditions,  S-matrix  elements  are  presented  for  reflection 
from  absorbing  boundary  conditions  for  the  light-light-light,  medium-light-medium 
and  heavy-light-heavy  mass  configurations.  Chapter  VI  presents  S-matrix  elements 
using  the  model  two  dimensional  system  of  two  coupled  Morse  oscillators  for  a  light- 
heavy-light  mass  configuration.  S-matrix  elements  for  this  mass  configuration  do 
not  exhibit  errors  attributed  to  reflection  from  absorbing  boundary  conditions.  An 
analysis  of  the  effect  of  the  depth  of  the  well  in  the  interaction  region  is  presented 
for  the  light-heavy-light  mass  configuration  as  well  as  the  effects  of  kinetically  de¬ 
coupling  the  two  Morse  oscillators.  Conclusions  and  recommendations  are  contained 
in  Chapter  VII.  Appendix  A  contains  the  derivation  of  the  channel  packet  method 
formulation  of  S-matrix  elements  for  two  dimensional  scattering.  Appendix  B  gives 
the  derivation  of  the  momentum  transformation  between  Jacobi  and  bond  repre¬ 
sentations.  Appendix  C  briefly  overviews  the  one  dimensional  and  two  dimensional 
codes  implementing  the  channel  packet  method  with  absorbing  boundary  conditions. 
Appendix  D  contains  the  complete  set  of  graphs  for  the  research  presented  in  Chap¬ 
ter  VI. 

1.4  Conventions  and  units 

Throughout  this  document,  the  term  channel  is  used  to  denote  a  given 
arrangement  of  atoms,  for  example,  A  -f  BC  and  AC  +  B.  The  term  channel  is  syn¬ 
onymous  with  arrangement  and  arrangement  channel.  Within  a  given  channel,  the 
fragments  involved  in  reactive  scattering  may  have  internal  degrees  of  freedom.  In 
quantum  mechanics,  these  internal  degrees  of  freedom  are  described  by  eigenstates 
and  associated  quantum  numbers.  The  term  subchannel  is  used  to  denote  the  com¬ 
bination  of  a  channel  with  a  particular  set  of  eigenstates  of  the  internal  degrees  of 
freedom  in  that  channel. 
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The  units  used  in  the  actual  computations  are  atomic  units  denoted  by  au. 
Table  1.1  contains  conversions  between  au  and  other  commonly  used  units. 

Table  1.1  Atomic  units 


unit 

1  au  equals 

mass 

9.11x10"^^  kg 

angular  momentum  {h) 

1.05x10“^'*  J-sec 

length 

5.29x10-11  m 

time 

2.42x10-1^  sec 

energy 

27.21  eV 
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II.  Quantum  scattering  theory 

The  central  goal  of  quantum  scattering  theory  is  to  determine  the  probabil¬ 
ity  of  reaction.  In  this  chapter,  the  theories  of  coordinate  systems,  Hamiltonians  and 
asymptotic  limits,  absorbing  boundary  conditions,  and  the  channel  packet  method 
and  the  application  of  Mpller  Operators  are  developed. 

2.1  Background 

Quantum  scattering  begins  with  the  time  dependent  Schrodinger  equation, 

=  (2-') 

where  H  is  the  Hamiltonian  and  the  state  of  the  system  is  given  by  (Q).  For  time 
independent  Hamiltonians,  the  formal  solution  to  Equation  2.1  is 

(t))  =  U  (t)  {t  =  0))  =  exp 

where  U  (t)  is  the  tirne  evolution  operator.  The  Hamiltonian  H  is  usually  defined  as 
the  sum  of  two  operators, 

H  =  f  +  V,  (2.3) 

where  f  is  the  kinetic  energy  operator  and  V  is  the  potential  energy  operator.  For 
a  wide  variety  of  reactions,  f  is  the  same  and  the  dynamics  are  governed  primarily 
by  y. 

2.2  Coordinate  systems 

Equation  2.1  must  be  expressed  in  some  representation  in  order  to  be  solved 
and  two  coordinate  systems  used  to  provide  a  representation  are  of  interest:  Bond 
coordinates  and  Jacobi  coordinates.'‘2  Section  2.3,  we  will  explore  the  need  for  two 
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coordinate  systems  as  they  relate  to  the  Hamiltonian  and  asymptotic  limits.  Bond 
coordinates,  illustrated  in  Figure  2.1,  treat  each  channel  on  an  equivalent  basis.  Two 
sets  of  Jacobi  coordinates,  {vo,,  Ra)  and  {rj^^R/s),  illustrated  in  Figure  2.2,  are  used 
to  provide  a  separable  representation  of  H  in  the  asymptotic  limit  in  each  of  the  two 
possible  channels  A  +  BC  and  AB  +  C . 

In  the  A  +  BC  channel,  the  transformation  between  Jacobi  coordinates  and 
bond  coordinates  in  two  dimensions  is  given  by 


where  rtij  is  the  mass  of  the  atom  j.  From  Equation  2.4,  Figure  2.1  and  Figure  2.2, 
it  is  clear  that  the  Jacobi  coordinate  Ra  is  equal  to  the  bond  coordinate  Y  and  ra 
is  equal  to  the  sum  of  and  distance  from  atom  B  to  the  center-of-mass  of  the  BC 
diatom. 

2.3  Hamiltonians  and  asymptotic  limits 

The  issue  of  which  coordinate  system  to  use  is  closely  linked  to  the  issue 
of  asymptotic  limits.  In  quantum  scattering,  we  must  assume  that  at  some  point  in 
space  the  fragments,  either  products  or  reactants,  are  widely  separated  such  that  the 
interaction  potential  between  them  is  zero.  Under  this  assumption,  the  Hamiltonian 
is  written  as 

H  =  Ho +  V,  (2.5) 

where  H  is  the  full  scattering  Hamiltonian,  Hq  is  the  asymptotic  Hamiltonian  and 
V  is  the  interaction  potential  that  vanishes  in  the  asymptotic  limit, 

lim  y  =  0,  7  =  a,^,  (2.6) 


where  7  =  a,  /3  labels  the  channel. 
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Figure  2.2  Jacobi  coordinates. 
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Using  Jacobi  coordinates,  the  asymptotic  Hamiltonian  can  be  separated  into 
parts  according  to 

Ho  =  Hrel  +  H,nt,  (2.7) 

where  Hrei  is  Hamiltonian  governing  the  relative  motion  between  reactants  or  prod¬ 
ucts  and  Hint  is  the  Hamiltonian  governing  the  internal  dynamics  of  reactants  or 
products.  It  is  important  to  note  that  Jacobi  coordinates  suitable  for  one  channel 
will  not  separate  Hq  in  another  channel.  .A.n  appropriate  choice  of  Jacobi  coordinates 
is  used  to  construct  the  initial  reactant  and  product  wave  packet  in  each  channel.  Ja¬ 
cobi  coordinates  are  used  since  the  Hamiltonian  is  separable  in  Jacobi  coordinates  as 
shown  in  Equation  2.7.  The  wave  packets  are  then  transformed  to  bond  coordinates 
to  facilitate  the  computation  of  a  correlation  function  used  to  calculate  S-matrix 
elements. 

2.4  Absorbing  boundary  conditions 

One  problem  with  time  dependent  methods  is  reflection  of  the  wave  function 
at  the  edges  of  the  grid  used  in  numerical  calculations.  These  reflections  can  lead 
to  spurious  results  and  are  usually  handled  by  choosing  a  grid  large  enough  that,  by 
the  time  the  w'ave  function  reaches  the  edge,  it  has  completely  left  the  interaction 
region  and  the  computation  may  be  terminated.  The  increased  grid  size  translates 
into  increased  computational  overhead.  Absorbing  boundary  conditions  have  been 
used  previously  in  a  variety  of  calculations^^"^^’^‘’^^’^®’^®’'‘°’^^^^  to  prevent  evolving 
wave  packets  from  reflecting  from  the  edge  of  the  grid  before  the  computation  has 
finished.  It  should  be  noted,  however,  that  reflection  can  only  be  minimized  and  not 
eliminated.  In  the  channel  packet  method,  it  is  in  the  calculation  of  the  correlation 
function  that  the  application  of  absorbing  boundary  conditions  proves  useful.  This 
application  is  discussed  in  Section  2.6.2. 

Negative,  imaginary  potentials  are  well  known  as  the  only  absorbing  potentials. 
These  potentials,  whose  mathematical  form  is  arbitrary,  absorb  by  exponentially 


2-4 


attenuating  the  wave  function  as  it  evolves  in  time.  In  the  application  of  absorbing 
boundary  conditions  to  the  channel  packet  method  performed  for  this  study,  a  simple 
exponential  of  the  form 

Va  -  ±.4iexp|^^^-^^| ,  (2.8) 

was  selected  where  A  and  B  were  chosen  to  minimize  reflection.  This  form  for  the 
absorbing  boundary  condition  was  chosen  because  it  is  a  smooth  function  which  can 
reach  a  large,  positive  value  in  a  relatively  short  distance.  These  two  characteris¬ 
tics  permit  the  use  of  the  smallest  grid  possible  with  a  subsequent  improvement  in 
computational  efficiency. 

2.5  Channel  packet  method 

The  time  dependent  channel  packet  method,  referred  to  as  the  channel 
packet  method,  was  developed  by  Weeks  and  Tannor^'^’^^’^*  and  is  based  on  the  use 
of  Mpller  operators.®^  In  this  section,  Mpller  operators  are  discussed  and  then  applied 
in  the  derivation  of  the  channel  packet  method. 

2.5.1  Mpller  operators.  The  development  of  Mpller  operators  begins 

with  the  Hamiltonian  and  asymptotic  limits  discussed  in  Section  2.3.  Utilizing  these 
ideas,  the  Hamiltonian  can  be  written  as 


H  =  +  (2.9) 

where  H  is  the  full  Hamiltonian,  Hq  is  the  asymptotic  Hamiltonian  in  channel  7 
and  V'^  is  the  interaction  potential  that  vanishes  in  the  asymptotic  limit.  The  label 
7  now  includes  the  channel  label  as  well  as  the  labels  of  the  internal  eigenstates, 
thereby  labeling  both  the  channel  and  subchannel. 
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M0ller  operators  are  defined  by 


fil  =  Jm  exp  (^.fj  exp  ’  l2-“» 

where  U  (t)  =  exp  ^—2^^  is  the  time  evolution  operator.  Mpller  operators  are 
isometric 

nlfii  =  1.  (2.11) 

and  if  bound  states  are  not  present,  they  are  unitary®^ 

=  (2.12) 

The  result  of  propagating  a  state  under  a  Mpller  operator  is  denoted  by 

^±i'i'in(out)}  =  1^±),  (2-13) 

where  |^±)  is  referred  to  as  a  Mpller  state.  The  effect  of  the  Mpller  operator  Q± 
is  to  propagate  the  state  |’^)  backward  (forward)  to  time  t  =  — oo  (t  =  Too)  under 
the  influence  of  ffo  ,  then  propagate  the  result  forward  (backward)  to  time  t  =  0 
under  the  influence  of  H.  By  definition,  the  Moller  operator  Cl-  operates  only 
on  asymptotic  product  states  l^ouf)  and  the  Mpller  operator  Cl+  operates  only  on 
asymptotic  reactant  states  Mpller  operators  are  used  to  define  the  scattering 

operator  through 

S  =  Cl[Cl+,  (2.14) 

where  the  probability  of  scattering  from  an  initial  state  l^m)  to  a  final  state  \’^out) 
is  given  by 

P  =  =  l{»-l«+>p.  (2.15) 
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Equation  2.15  yields  S-matrix  elements  in  the  out)  representation.  While 

this  is  an  acceptable  result,  it  is  more  useful  to  express  S-matrix  elements  in  the 
energy  or  momentum  representation.  In  the  channel  packet  method,  expressing  S- 
matrix  elements  in  the  momentum  representation  leads  to  the  calculation  of  a  time 
dependent  correlation  function. 

2.5.2  The  channel  packet  method.  The  core  of  the  channel  packet 

method  is  the  use  of  two  wave  packets  to  compute  state-to-state  S-matrix  elements. 
One  wave  packet  represents  the  asymptotic  reactant  subchannel,  labeled  7,  and  the 
second  wave  packet  represents  the  asymptotic  product  subchannel,  labeled  7'.  The 
choice  of  the  two  subchannels,  reactant  and  product,  determines  which  S-matrix 
elements  are  calculated.  The  technique  is  to  propagate  each  wave  packet  using  the 
appropriate  Mpller  operator. 

The  7  subchannel  asymptotic  Hilbert  space  is  spanned  by  the  direct  product 
of  eigenstates  of  the  relative  Hamiltonian,  represented  by  the  ket  |A:^),  and 

eigenstates  of  the  internal  Hamiltonian,  Hi^t,  represented  by  the  ket  I7).  The  basis 
vectors  spanning  the  Hilbert  space  are  given  by, 

\k„i).  (2.16) 

These  basis  vectors  are  eigenstates  of  the  asymptotic  Hamiltonian  with  eigenvalues 
given  by 

H^\k„'r)  =  {Hh  +  flli)  \Kl)  =  =  ^1^7,7),  (2-17) 

where  h^k'^j2pL..^  is  the  relative  kinetic  energy,  p-f  is  the  reduced  mass  of  the  system, 
is  the  energy  associated  with  the  internal  eigenstate  and  E  is  the  total  energy. 
A  new  set  of  states  is  generated  by  operating  on  the  basis  \k.y,'y)  with  the  Mpller 
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operators 


(2.18) 


where  the  ±  in  lA;-^,7±)  labels  both  the  state  and  the  state  I7). 

It  is  easy  to  show  that  the  l/L-y,7±)  are  eigenstates  of  the  full  Hamiltonian. 
Using  the  intertwining  relationship,^^ 

M±=U±Ho,  (2.19) 

and  from  Equations  2.17  and  2.18,  it  follows  that 

H\K^±)  =  HQl\k,,^) 

=  +  (2.20) 

where  the  eigenvalues  of  H  are  the  same  as  the  eigenvalues  of  Hq. 

Consider  the  initial  state  |^7n(out))  constructed  as  a  product  of  a  linear  combi¬ 
nation  of  kets  [A;-,.)  with  an  eigenstate  I7), 

/+00 

dk,,li(k,)\Kt)-  (2.21) 

“OO 

The  corresponding  Mpller  states  |^^±)  are  given  by 

/+00 

rfM±(fe,)|i,.7±>.  (2.22) 
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From  Equations  2.21  and  2.22  it  is  clear  that  the  expansion  coefficients, 
used  to  expand  the  initial  wave  function  in  terms  of  the  eigenstates  of  Hq  are  the 
same  coefficients  used  to  expand  the  Mpller  states  in  terms  of  the  eigenstates  of 
H.  Thus,  the  eigenstates  of  the  full  Hamiltonian  can  be  used  to  expand  the  Mpller 
states,  |^±). 

The  isometric  subchannel  Mpller  operators,  Oj.,  acting  on  the  basis  vectors 
Ik'Yi'y),  which  span  the  asymptotic  Hilbert  space,  generate  two  sets  of  basis  vectors 
|A:-^,7±),  which  are  eigenvectors  of  the  full  scattering  Hamiltonian  H.  The  two 
complete  sets  satisfy  the  orthogonality  relations,®^ 


and 


=  {kl,,7'\k^,l) 


(2.23) 


S{E'-E)Sl\ 


(2.24) 


where  ^\k^\  —  is  the  density  of  states  and  is  the  on-shell  S-matrix  ele¬ 

ment  giving  the  probability  amplitude  of  scattering  from  subchannel  7  with  relative 
momentum  k^  to  subchannel  7'  with  relative  momentum  fcy  . 

Equation  2.24  shows  that  all  that  is  necessary  to  obtain  S-matrix  elements  in 
the  momentum  representation  is  the  evaluation  of  the  scalar  product  {ky,Y—\k-y,  7-t-). 
In  the  channel  packet  method  this  is  done  by  first  considering  the  Fourier  transform 
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of  the  time  evolution  of  the  Moller  state  | 


Substituting  the  expansion  of  |^+)  in  terms  of  |A:^,7),  given  in  Equation  2.22,  into 
Equation  2.25,  it  is  easy  to  show  that  lA].  {£))  is  an  unnormalized  eigenvector,  shown 
in  Appendix  A,  of  the  full  scattering  Hamiltonian  given  by 

lAX  (£))  =  (+*:,)  l+fc,, 7+)  +  <1-  (-k,)  I-*,, 7+)).  (2.26) 

The  prefix  +  or  -  on  fc^  in  Equation  2.26  explicitly  labels  the  degenerate  eigenstates 
of  the  full  scattering  Hamiltonian.  The  degenerate  eigenstates  are  explicitly  labeled 
by  the  ±  prefix  as  having  positive  momentum  (+)  or  negative  momentum  (  — ).  S- 
matrix  elements  result  from  the  evaluation  of  the  scalar  product  (^1  \A\{E)),  where 
Equation  2.22  is  used  to  expand  I'I'l’)  in  terms  of 

+ll  (k'y)  7+  (-«:,)  (y.,7'-|-<:„7+)|.  (2.27) 

The  orthogonality  relation  given  in  Equation  2.24  reduces  the  integral  in  Equa¬ 
tion  2.27  to 


+  rjZ  (+K')  ‘^+v,-fc7 

+  ri‘_  (-fc;,)  7+  i+k-,)  5!^ 

+  (~^v)  *7+  ("^7)  SZk'^_k^]- 


(2.28) 


We  are  free  to  simplify  Equation  2.28  since  the  only  requirement  that  the  asymptotic 
states  |^7n(ouq)  na^st  satisfy  in  order  to  obtain  Equation  2.28  is  that  they  are  elements 


2-10 


of  the  7*^  subchannel  Hilbert  space.  The  simplification  is  to  choose  the  expansion 
coefficients  r)±  {k-y)  such  that  they  include  only  one  degenerate  state  for  each  value 
of  the  relative  kinetic  energy  in  the  7  subchannel.  Choosing  7+  {k~^)  =  0  for  >  0 
and  7_  =  0  for  k'^  <  0,  the  evaluation  of  Equation  2.28  reduces  to 


=  til  {+K')  '?+(■ 


(2.29) 


Now,  we  make  a  slight  notation  change  where  the  label  ±  on  the  relative 
momentum  corresponds  to  the  four  possible  combinations  of  |  ±  'f7n(ou<))-  Inverting 
Equation  2.29,  we  obtain  the  on-shell  S-matrix  element  for  scattering  from  the  state 
\  —  k^^'f)  to  the  state  |  +  k-y',Y), 


07  W 
7 


(27r)  ^  h 

n-  {+K')  ’i* 


(2.30) 


where  the  scalar  product  (E))  is  given  by  the  Eourier  transform  of  the  cor¬ 
relation  function  (^1  jexp  1^+); 


{r:\AiiE))  = 


r-H-co  ,  fiEt\  ,.,v, 
dt  exp  \Y~^J 


in)- 


(2.31) 


The  numerical  implementation  of  Equations  2.30  and  2.31  will  be  covered  as  part  of 
the  next  section. 
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2.6  Numerical  implementation 

The  numerical  implementation  of  the  channel  packet  method  is  fairly 
straight  forward.  In  this  section,  we  will  outline  the  steps  involved  in  computing 
S-matrix  elements  for  the  collinear  reaction 

T  +  BC  (u)  ^  AB  {N)  +  C  (2.32) 

where  the  internal  diatom  eigenstates  are  labeled  v  and  v'.  The  actual  calculation 
of  S-matrix  elements  can  be  broken  into  three  broad  sections:  calculating  the  Mpller 
states,  calculating  the  correlation  function  and  calculating  S-matrix  elements. 

2.6.1  Calculating  the  Mpller  states.  Numerically,  the  reactant  Mpller 

state  is  calculated  in  the  following  manner.  In  the  asymptotic  limit  of  subchannel  7, 
the  asymptotic  channel  Hamiltonian  is  separable  using  Jacobi  coordinates  according 
to  Equation  2.7.  is  a  function  only  of  Ry  and  is  a  function  only  of  Vy.  The 
solutions  to  are  diatomic  vibrational  eigenstates  and  the  solutions  to 
plane  waves.  The  reactant  state  |^7n)  is  the  direct  product  of  a  single  vibrational 
eigenstate  and  a  linear  combination  of  eigenstates  of  given  by  Equation  2.21. 
The  reactant  state  j^/n)  is  then  transformed  to  bond  coordinates  for  propagation. 
The  choice  of  t  =  0  and  the  placement  of  the  wave  functions  is  purely  arbitrary, 
though,  the  choice  of  ^  =  0  and  placement  in  the  interaction  region  facilitate  the 
calculation  of  the  correlation  function  in  Section  2.6.2. 

The  state  is  propagated  using  the  Mpller  operator  as  follows.  The 

limits  t  ±00  are  numerically  approximated  by  t  ->■  ±r  where  the  time  r  is  taken 
to  be  when  the  wave  packet  completely  exits  the  interaction  region.  The  interaction 
region  is  taken  to  be  where  is  not  negligibly  small.  First,  |’^7n)  is  propagated 
backward  in  time  from  t  =  0  to  t  =  — r  under  Hq.  This  can  be  done  analytically 
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since  the  internal  eigenstate  I7)  evolves  in  time  as 


I7  (0)  =  exp  I7  {t  =  0))> 


(2.33) 


and  the  time  evolution  of  ri{k)  is  given  by, 


53 


ri{k)  =  [  dx  (2t:  [Axo)^] 

J  —  00  \  I 


iht 


2fi  (Axo) 


X  exp  < 


0  2m 


1  + 


iht 


2(AxorM 


exT[>{ikx).  (2.34) 


Equation  2.34  is  a  complex  Gaussian  where  Axq  is  the  initial  wave  packet  spread  in 
the  coordinate  representation  and  ko  is  the  initial  momentum  offset.  The  advantage 
to  using  Equations  2.33  and  2.34  is  that  half  the  action  of  the  Moller  operator  fi±, 
the  evolution  of  the  initial  state  under  Hq,  can  be  done  analytically. 

The  numerical  propagation  is  based  on  the  Baker-Campbell-Housdorf  theo- 
rem^^  which  states 


exp  (^A+  -  exp  exp  exp  B])  (2.35) 

for  operators  A  and  B.  For  the  time  evolution  operator  U  (f),  applying  the  theorem 
yields 

exp[^  (f  +  V)Af]  =  exp  Qf  Af^  exp  exp  ^-[f,V]At^'j  ,  (2.36) 

where  f  is  the  kinetic  energy  operator  and  V  is  the  potential  energy  operator.  Simply 
neglecting  the  last  term  leads  to  an  error  on  the  order  of  Af^  .  Instead,  we  use  the 
split  operator  approach  to  achieve  the  same  effect  with  improved  accuracy.^®  The 
split  operator  method  first  splits  the  Hamiltonian  in  half  before  expanding  in  terms 
of  f  and  V.  Then  the  Baker-Campbell-Housdorf  is  applied  to  the  two  exponentials. 
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This  leads  to  an  approximation  of  Equation  2.36  given  by 

exp[^  (f  +  «  exp  exp  ^^TAtj  e.xp  +  O  (At^)  , 

(2.37) 

where  the  error  in  the  approximation  is  on  the  order  of  At^  instead  of  At^.  In 
Equation  2.37,  the  potential  operator  is  diagonal  in  the  coordinate  representation 
and  the  kinetic  energy  operator  is  diagonal  in  the  momentum  representation.  The 
transformation  between  the  two  representations  is  a  Fourier  transform  pair  given  by 

g{k)  =  J^~^{g{x)}  = -^  [  dxexp{i2Trkx)g{x),  (2.38) 

v27r  J-oo 

and 

1  /'+00 

g{x)  =  J^{g{k)}  = -^  dkexp{-i2TTkx)g{k).  (2.39) 

Numerical  propagation  proceeds  by  first  defining  a  coordinate  grid  and  a  momen¬ 
tum  grid  using  the  bond  coordinate  representation.  The  coordinate  grid  is  used  to 
determine  the  diagonal  representation  of  exp  At^  and  the  momentum  grid  is 
used  to  determine  the  diagonal  representation  of  exp  (^^TAt^.  The  initial  reactant 
wave  packet  '5  (r-^,  Ity.t  =  0)  is  transformed  to  the  bond  coordinate  representation 
and  using  Equation  2.40  the  wave  packet  at  f  -t-  Af  is  computed, 

1^1  ((  +  AO)  =  exp  (^VAf)  :r{exp  (ifAi)  :r-'{exp  (^V-Ai)  (0)}}. 

(2.40) 

In  order  to  efficiently  propagate  using  Equation  2.40,  we  use  Fast  Fourier  transforms 
(FFTs).  Repeated  application  of  Equation  2.40  iteratively  determines  the  wave 
packet  at  time  r.  The  calculation  of  the  product  state  is  similar. 

2.6.2  Calculating  the  correlation  function.  With  the  two  Mpller  states 

and  !'!’_)  calculated  at  f  =  0,  the  next  step  is  to  compute  the  correlation 


2-14 


function  in  Equation  2.31  and  given  by 

Cv7(0  =  (^-|e-^p  1^+)- 

We  need  to  know  Cy-^  (0  from  time  <  =  — r  to  t  =  r  in  order  to  take  the  Fourier  trans¬ 
form  in  Equation  2.31.  First,  the  reactant  Mqller  state  l’F+)  is  propagated  from  time 
f  z=  0  to  t  =  r  under  the  full  scattering  Hamiltonian  H.  Second,  the  reactant  M0ller 
state  !'!'+)  at  t  =  0,  is  propagated  from  time  t  =  0  to  t  =  — r  under  the  full  scattering 
Hamiltonian  H  again  computing  the  correlation  at  each  step.  It  is  in  calculating  the 
correlation  function  that  absorbing  boundary  conditions  are  applied  since  the  Mqller 
states  are  typically  well  localized  in  the  interaction  region  at  t  =  0.  The  absorbing 
boundary  conditions  are  applied  in  an  area  where  {t  =  0))  is  already  zero  thus 
C-f'-y  is  zero  and  the  absorbing  boundary  conditions  do  not  affect  its  value.  The  form 
of  the  absorbing  boundary  conditions,  on  the  discretized  numeric  grid,  are  altered 
slightly  from  those  in  Equation  2.8  by  the  multiplication  of  a  step  function  so  that 
the  absorbing  boundary  conditions  are  in  fact  zero  where  |'Fl  [t  =  0))  is  non-zero. 
The  absorbing  boundary  conditions,  as  implemented  numerically  are  given  by 

Va  -  ±h{x  -  aro)Azexp  ,  (2.42) 

where  h{x  —  xq)  is  the  heavyside  step  function  given  by 

^  1.  Xn  X  ^ 

h{x-xo)=  .  •  (2-43) 

0,  otherwise  j 

2.6.3  Calculating  S-matrix  elements.  In  the  channel  packet  method, 

the  formula  for  computing  S-matrix  elements  is  contained  in  Equations  2.30  and 


2-15 


2.31.  Numerically,  the  starting  point  is  Equation  2.31  given  by 


/+OC  /iEt\ 

dtexpy—jCy^it),  (2.44) 


where  C~f'-y{t)  is  the  correlation  function  computed  in  the  previous  section.  Numeri¬ 
cally,  Equation  2.44  is  implemented  by  taking  the  discrete  Fourier  transform  of  the 
correlation  function  computed  in  Section  2.6.2.  We  are  free  to  chose  the  grid  on  which 
the  discrete  Fourier  transform  is  taken.  The  best  choice  is  an  energy  grid  where  the 
underlying  momenta,  recall  E  =  are  such  that  the  expansion  coefficients  r]±  in 
Equation  2.21  are  numerically  non-zero.  This  ensures  that  the  normalization  by  the 
expansion  coefficients  in  the  numeric  implementation  of  Equation  2.30  does  not  lead 
to  numeric  error  caused  by  dividing  large  numbers  by  small  numbers.  After  taking 
the  discrete  Fourier  transform  of  the  correlation  function,  the  computation  of  S- 
matrix  elements  using  Equation  2.30  is  simple  and  straight  forward.  The  probability 
of  reaction  is  given  by 


p^'y 

7 


(2.45) 


where  are  the  S-matrix  elements  computed  from  Equation  2.30. 

7' 


2-16 


III.  One  dimensional  quantum  scattering 

As  an  introduction  to  quantum  scattering,  the  channel  packet  method  with 
absorbing  boundary  conditions  was  applied  to  quantum  scattering  in  one  dimension 
for  two  systems:  A  square  well,  and  a  Gaussian  well  with  Gaussian  barriers. 


3.1  One  dimensional  square  well 

The  computation  of  S-matrix  elements  begins  with  Schrodinger’s  Equation 
in  one  dimension. 


(3.1) 


where  the  Hamiltonian  H  is  given  by 


H  =  T  +  V, 


(3.2) 


f  is  the  kinetic  energy  operator  and  V  is  the  potential  operator.  Since  there  is  only 
a  single  degree  of  freedom,  the  particle  has  no  internal  structure.  The  formula  for  S- 
matrix  elements  using  the  channel  packet  method,  given  in  Equation  2.30,  simplifies 
to 

^/\k'\\k\ 

(3.3) 


^  (27r)  'ft  I, 


where  the  scalar  product  (^_|A+  (E))  is  given  by  the  Fourier  transform  of  the  cor¬ 
relation  function 

/•H-oo  /  Et\  (  Ht\ 

{».|At(£))=  *exp[iyj(«.|exp(-.— (3.4) 

The  mathematical  form  of  the  one  dimensional  square  well  is  given  by 


V{x)  = 


—Vo,  —a  <  X  <  a 
0,  otherwise 


(3.5) 
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where  the  potential  V'(a;)  is  non-zero  from  a;  =  — a  to  x  =  a.  This  square  well 
is  shown  in  Figure  3.1.  It  should  be  noted  that  the  sides  of  the  one  dimensional 
square  w'ell  are  discontinuous  at  x  =  ia.  This  discontinuity  and  the  impact  on  the 
numerical  implementation  is  discussed  further  in  Section  3.1.1 

For  the  potential  illustrated  in  Figure  3.1,  the  scattering  is  elastic  and  S-matrix 
elements  are  non-zero  only  for  =  |A:|.  For  a  one  dimensional  square  well  potential 
energy  surface  the  Hamiltonian  is  given  by 


- 


—Vq,  —a<x<a 
0,  otherwise 


(3.6) 


and  the  probability  of  transmission,  given  by  Pk>k  =  has  an  analytic  solu¬ 

tion^®’’"  given  by 


P{E)  =  [cos^  {2ka)  + 


(3.7) 


with 

k  = 
k'  = 

and  d  = 

where  2a  is  the  well  width,  Vo  is  the  we 
particle  mass. 


/|?(£-V4) 

(3.8) 

i/lS 

(3.9) 

k  k' 
k''^  k' 

(3.10) 

depth,  E  is  the  particle  energy  and  m  the 


3.1.1  Choosing  the  grid.  The  first  step  in  applying  the  channel  packet 

method  to  a  given  scattering  problem  is  determining  the  grid  parameters.  Kosloff’s 
seminal  article’®  on  time  dependent  scattering  outlines  the  constraints  on  choosing 
grids.  Consider  first  the  total  energy  a  particle  may  have  during  the  scattering  event. 
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Figure  3.1  One  dimensional  square  well  of  width  2a  and  depth  Vq. 
The  total  energy,  Fimax,  given  by 


E 


max 


Eint  +  Tlnax  “  Kr 


(3.11) 


is  not  the  same  as  the  expectation  value  of  the  Hamiltonian  (H)-  Instead,  it  is  a 
combination  of  the  maximum  translational  energy,  the  maximum  potential  energy 
and  the  internal  energy.  In  one  dimensional  scattering,  there  is  no  internal  energy 
to  consider,  leaving  only  translational  and  potential  energy.  The  maximum  transla¬ 
tional  energy  is  taken  to  be  the  energy  of  the  single  component  of  the  wave  packet 
with  the  highest  absolute  value  of  k.  For  the  potential  energy,  the  maximum  depth 
of  the  potential  well  is  added  to  the  total  energy  in  Equation  3.11.  Barrier  energies 
are  not  subtracted  from  the  total  energy. 

From  the  total  energy,  the  maximum  value  of  the  momentum  grid,  fcmax  is 
derived  through  the  relationship  between  energy  and  momentum  for  non-relativistic 
scattering, 

fi^max  =  \/2mEmax-  (3.12) 
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Figure  3.1  One  dimensional  square  well  of  width  2a  and  depth  Vb- 
The  total  energy,  £^max,  given  by 


Eint  +  Tma.x  ~  Kn 


(3.11) 


is  not  the  same  as  the  expectation  value  of  the  Hamiltonian  (H).  Instead,  it  is  a 
combination  of  the  maximum  translational  energy,  the  maximum  potential  energy 
and  the  internal  energy.  In  one  dimensional  scattering,  there  is  no  internal  energy 
to  consider,  leaving  only  translational  and  potential  energy.  The  maximum  transla¬ 
tional  energy  is  taken  to  be  the  energy  of  the  single  component  of  the  wave  packet 
with  the  highest  absolute  value  of  k.  For  the  potential  energy,  the  maximum  depth 
of  the  potential  well  is  added  to  the  total  energy  in  Equation  3.11.  Barrier  energies 
are  not  subtracted  from  the  total  energy. 

From  the  total  energy,  the  maximum  value  of  the  momentum  grid,  is 
derived  through  the  relationship  between  energy  and  momentum  for  non-relativist ic 
scattering, 

hkm^x  =  y/^rnE^.  (3.12) 
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The  momentum  and  coordinate  representations  are  related  through  a  Fourier  trans¬ 
form  pair,  given  by  Equations  2.38  and  2.39.  However,  this  choice  of  Fourier  trans¬ 
form  pair  requires  the  channel  packet  method  to  make  use  of  a  different  momentum 
usually  denoted  b\^  h  and  given  by 

h  =  (3.13) 

27r 

Throughout  this  document,  kb  is  used  and  the  subscript  b  is  omitted.  The  Fourier 
transform  pair  in  Equations  2.38  and  2.39  also  yields  the  relationship  between  /cmax 
and  Ax  the  coordinate  grid  spacing, 

Ax  =  -7^^ — .  (3.14) 

-^max 

The  choice  of  the  maximum  value  of  the  coordinate  grid  is  less  straightforward.  The 
usual  manner  in  which  x^ax  is  determined  is  to  consider  the  time  evolution  of  the 
Mpller  states.  The  coordinate  grid  must  be  large  enough  to  support  the  evolving 
Mpller  states  until  the  calculation  of  the  correlation  function  can  be  terminated. 
The  purpose  behind  using  absorbing  boundary  conditions  is  to  decrease  the  size  of 
the  grid  and  exponentially  attenuate  the  evolving  Mpller  state  before  it  reaches  the 
edges  of  the  smaller  grid.  Thus,  the  grid  must  only  be  large  enough  to  contain  the 
interaction  region  “  as  well  as  the  Mpller  states  at  t  =  0  and  the  absorbing  boundary 
conditions.  The  absorbing  boundary  conditions  are  placed  as  close  as  possible  to 
the  Mpller  states  without  overlapping  them.  This  maximizes  the  benefit  of  using 
absorbing  boundary  conditions.  The  extent  of  the  grid  past  the  point  at  which 
the  absorbing  boundary  conditions  are  placed  need  only  be  large  enough  so  that 
the  evolving  Mpller  state  is  fully  attenuated  before  it  reflects  from  the  edge  of  the 
grid.  Initially,  the  location  of  the  absorbing  boundary  conditions  is  based  upon  the 
estimated  spatial  extent  of  the  Mpller  states  with  the  absorbing  boundary  conditions 

“The  area  in  which  the  interaction  potential  is  non-zero. 
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location  adjusted  later  as  necessary.  The  numerical  implementation  of  the  channel 
packet  method,  discussed  in  Section  2.6,  makes  use  of  FFTs  and  the  Temperton'’’ 
FFT  used  in  this  project  requires  that  be  a  power  of  2.  3  or  5.  The  most  efficient 
FFT  is  for  N  to  be  a  power  of  2.  In  choosing  to  be  a  power  of  2,  that  restriction 
combined  with  the  estimated  x^ax  leads  to  choosing  N  such  that  it  is  a  power  of  2 
and  Y  x  Ax  >  Xmax-  In  one  dimension,  we  use  y  since  the  spatial  grid  contains  both 
positive  and  negative  values  of  x  and  is  chosen  to  be  symmetric. 

The  total  energy  in  Equation  3.11  is  also  used  to  determine  the  time  step  At. 
Similar  to  the  momentum  and  coordinate  representations,  energy  and  time  are  also 
related  by  a  Fourier  transform  pair  which  yields  the  relationship  between  Emax  and 

Atmax  =  - •  (3.15) 

^max 

At  this  point,  the  parameters  Atmax)  ^max?  and  Ax  have  been  determined. 

The  one  dimensional  square  well  is  discontinuous  at  x  =  ±a  as  shown  in 
Equation  3.5  and  Figure  3.1.  Using  a  discrete  spatial  variable  x  leads  to  a  modeling 
of  the  discontinuous  sides  as  slopes  instead  of  perfectly  vertical.  Numerically,  it  would 
require  an  infinite  number  of  points  in  the  spatial  grid  to  reach  vertical  sides.  Fewer 
points  leads  to  the  square  well  being  modeled  as  a  trapezoid.  As  the  number  of  points 
increases,  the  slope  of  the  sides  of  the  trapezoid  approaches  the  perfectly  vertical, 
and  discontinuous,  sides  of  the  square  well.  This  is  shown  in  Figure  3.2  where  the 
numeric  implementation  of  the  analytic  square  well  is  shown  as  the  number  of  grid 
points  increases  and  their  spacing  decreases.  With  an  infinite  number  of  points,  the 
numeric  potential  perfectly  matches  the  analytic  potential  in  Equation  3.5.  Using 
the  Ax  previously  determined,  the  number  of  points  N  needs  to  be  4096  in  order  to 
obtain  a  numeric  solution  that  is  close  to  the  analytic  solution. 

3.1.2  M0ller  states.  In  constructing  the  initial  wave  packets,  the 

reactant  channel  is  chosen  arbitrarily  to  be  the  left  side  of  the  potential,  x  <  0,  and 
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Figure  3.2  The  discontinuous  sides  of  the  one  dimensional  square  well  modeled  as 
trapezoidal  sides. 

the  product  channel  is  the  right  side  of  the  potential,  x  >  0.  Compact  Mpller  states 
are  important  in  reducing  the  overall  grid  size  and  the  most  compact  Mpller  states 
for  scattering  from  a  one  dimensional  square  well  are  those  with  an  initial  spread  in 
the  coordinate  representation  less  than  the  width  of  the  well.  This  determination 
of  compactness  is  based  solely  on  trial  and  error.  In  this  case,  the  reactant  Mpller 
state  was  computed  for  a  range  of  initial  spreads  and  then  plotted  on  the  spatial 
grid.  After  a  few  trials,  it  became  clear  that  the  narrower  spreads  were  leading  to 
more  compact  Mpller  states.  The  trials  were  then  performed  again  to  find  the  most 
compact  Mpller  state,  i.e.  narrowest  spatial  spread,  without  excessive  spread  in  the 
momentum  representation  as  the  spread  in  k  is  inversely  proportional  to  the  spread  in 
X.  The  two  wave  packets,  representing  reactants  and  products,  have  zero  coordinate 
offset  and  the  same  spread  in  the  coordinate  representation  in  order  to  locate  the 
compact  Mpller  states  in  the  interaction  region.  The  momentum  offset  must  be  large 
enough  so  that  the  initial  wave  packet,  in  the  momentum  representation,  has  only 
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positiv'e  or  nogcitive  momentum  in  order  to  apply  the  formula  for  S~matrix  elements 
given  in  Ecjuation  3.3.  With  the  initial  wave  packets  created,  the  reactant  M0llei 
state  is  computed  bv  propagating  the  initial  reactant  wave  packet  backward  in  time 
analytically  using  Equation  2.34  and  then  forward  in  time  numerically  under  the 
full  scattering  Hamiltonian  H  using  the  split  operator  approximation.  Likewise,  the 
product  Holier  state  is  computed  by  propagating  the  initial  product  wave  packet 
forward  in  time  analytically  and  then  backward  in  time  numerically  under  the  full 
scattering  Hamiltonian. 

3.1.3  Correlation  function.  In  calculating  the  correlation  function,  the 

reactant  Holler  state  remains  at  t  —  0  and  the  product  Holler  state  evolves  in  time. 
Eirst,  the  product  Holler  state  is  propagated  from  t  —  0  to  t  =  —t  while  computing 
the  correlation  function  at  every  step.  Second,  the  product  Holler  state  at  t  =  0, 
which  is  stored  in  memory,  is  propagated  to  t  =  +t,  again  computing  the  corre¬ 
lation  function  at  every  step.  Without  the  use  of  absorbing  boundary  conditions, 
the  coordinate  grid  would  have  to  be  large  enough  to  support  the  evolving  Holler 
state;  instead,  absorbing  boundary  conditions  are  used  to  exponentially  attenuate 
the  evolving  Holler  state  and  prevent  reflection  from  the  edges  of  the  smaller  grid. 
The  form  of  the  absorbing  boundary  conditions,  given  in  Equation  2.8,  needs  only 
three  parameters,  .4,  B ,  and  xq.  The  choice  of  xq  is  based  on  the  spatial  extent 
of  the  Holler  states  and  is  chosen  so  that  the  absorbing  boundary  conditions  be¬ 
gin  close  to  the  Holler  states  but  that  the  Holler  states  are  numerically  zero  at  xq- 
The  parameters  A  and  B  are  determined  by  the  need  to  fully  absorb  the  evolving 
Holler  state  while  minimizing  reflection.  Reflection  is  minimized  by  selecting  an 
absorbing  boundary  condition  amplitude  large  enough  to  absorb  the  highest  energy 
components  of  the  evolving  Holler  state.  It  is  also  important  to  avoid  the  sudden 
onset  of  absorbing  boundary  conditions  that  will  tend  to  reflect  the  lower  energy 
components  of  the  evolving  Holler  state.  As  investigated  in  Section  3.2.3,  absorbing 
boundary  conditions  that  reach  a  height  of  7  au  are  sufflcient  to  fully  attenuate  the 
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evolving  M0ller  state  before  it  reaches  the  edge  of  the  grid.  The  initial  step  has  only 
a  small  effect  and  should  be  in  the  range  of  1  x  10“^  to  1  x  10“^.  The  parameter  A 
is  determined  by  two  competing  factors:  the  need  to  place  the  absorbing  boundary 
conditions  in  as  small  a  space  as  possible  to  reduce  the  grid  size  while  avoiding  the 
steep  onset  of  absorbing  boundary  conditions.  A  safe  first  choice  for  A  is  on  the 
order  of  1  x  10“®.  While  a  larger  value  of  A  may  lead  to  more  compact  absorbing 
boundary  conditions,  it  may  possibly  lead  to  reflection  errors.  A  small  value  will 
most  likely  eliminate  that  possibility. 

3.1.4  Results  and  convergence  testing.  The  initial  parameters  of  the  one 

dimensional  square  well  and  the  resulting  grid  parameters,  are  listed  in  Table  3.1. 
The  small  time  step  and  the  fact  that  potential  wells  tend  to  trap  long-lived,  quasi¬ 
bound  states  requires  a  long  computation  time  in  order  for  the  correlation  function  to 
reach  zero  numerically.  The  choice  of  zero  is  important.  Choosing  to  terminate  too 
early  leads  to  incomplete  information  for  computing  S-matrix  elements  and  a  resul¬ 
tant  oscillation  in  the  probability  of  transmission.  Reviewing  the  literature  reveals 
that  a  typical  choice  for  numeric  zero  is  on  the  order  of  \C\  = 

This  is  the  cutoff  used  throughout  this  research  project.  The  reactant  Mpller  state 
is  shown  in  Figure  3.3  and  the  resulting  correlation  function  is  shown  in  Figure  3.4. 
The  probability  of  transmission,  both  numeric  and  analytic,  are  shown  in  Figure  3.5. 
They  are  in  agreement  and  show  that  the  combination  of  the  channel  packet  method 
with  absorbing  boundary  conditions  provides  a  valid  approach  to  computing  S- 
matrix  elements.  The  main  source  of  error  is  that  the  discontinuous  square  sides  of 
the  one  dimensional  square  well  are  computationally  modeled  as  trapezoids. 

The  performance  of  a  numerical  solution  needs  to  be  evaluated  against  two 
tests  of  convergence.  First,  does  the  solution  converge  to  the  correct  solution?  A 
solution  which  converges  to  a  wrong  answer  is  of  little  use.  The  channel  packet 
method  uses  discrete  spatial  and  temporal  variables  and  the  solution  is  only  approx¬ 
imate.  However,  as  the  spatial  and  temporal  grids  are  refined,  the  channel  packet 
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Table  3.1  Initial  parameters  for  scattering  from  a  one  dimensional  square  well. 


Parameter 

Value  (au) 

Notes 

mass 

1.0 

Electron  mass 

Ax 

0.005 

X  grid  resolution 

nstep 

450 

Number  of  time  steps  for  reactant  propagation 

npstep 

1000 

Number  of  time  steps  for  correlation  function  propagation 

At 

0.001 

Time  step 

ioff 

0.0  1 

initial  x  offset 

imom 

1.8 

initial  momentum  offset  in  units  of 

sprd 

0.25 

initial  wave  packet  spread  in  coordinate  representation 

X(au) 


Figure  3.3  Reactant  Mpller  state  for  a  one  dimensional  square  well.  The  well  ex¬ 
tends  from  —1  to  1. 
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Figure  3.4  Absolute  value  of  the  correlation  function  for  scattering  from  a  one 
dimensional  square  well. 


Figure  3.4  Probability  of  transmission  for  a  one  dimensional  square  well. 
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Figure  3.6 


Convergence  of  the  channel  packet  method  solution  towards  the  analytic 
solution  for  a  one  dimensional  square  well. 


method  should  converge  to  the  exact  solution.  The  Li  norm  is  used  as  a  measure  of 
convergence  and  is  defined  by 


N 


N 


E 


\Pj, analytic  Pj\ 
Pj, analytic 


(3.16) 


where  N  is  the  number  of  points  computed  in  the  probability  of  reaction  curve, 
Pj, analytic  is  the  analytic  solution  at  point  j  and  P,  is  the  calculated  probability.  The 
solution  is  considered  to  have  converged  when  the  Li  norm  has  converged  to  three 
digits.  For  the  one  dimensional  square  well,  this  is  illustrated  in  Figure  3.6  where 
the  spatial  grid  resolution  Ax  is  gradually  reduced  and  the  Li  norm  has  converged 
to  three  digits  and  the  method  is  considered  to  have  converged  to  the  correct  answer. 

The  second  test  of  convergence  illustrates  the  efficiency  of  the  channel  packet 
method.  This  test  concerns  how  the  approximate  solution  approaches  the  analytic 
solution  w'hen  one  grid  resolution  is  varied  while  the  rest  remain  fixed.  In  the  first 
case,  the  temporal  resolution  remains  fixed  while  the  spatial  resolution  is  changed. 
In  the  second  case,  the  temporal  resolution  is  varied  while  the  spatial  resolution 
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remains  fixed.  However,  there  is  one  proviso  to  the  case  where  the  spatial  resolution 
is  changed.  Since  the  one  dimensional  square  well,  in  the  approximate  solution, 
consists  of  trapezoidal  sides  ^  it  is  important  to  model  the  trapezoidal  sides  so  that 
the  improvement  in  the  accuracy  of  the  solution  is  due  solely  to  the  changes  in  the 
spatial  grid  and  not  the  trapezoidal  potential  approaching  the  analytic  square  well 
potential.  The  trapezoidal  result  of  discretizing  the  spatial  grid  was  discussed  in 
Section  3.1.1  and  illustrated  in  Figure  3.2.  In  determining  the  order  of  convergence 
for  changes  in  Ax,  the  one  dimensional  square  well  potential  given  in  Equation  3.5 
is  modified  to  explicitly  model  the  trapezoidal  potential.  This  new  potential  is 
illustrated  in  Figure  3.7  and  given  by 


V(x)  = 


0,  X  <  —do 

m2  X  X  +  62,  —do  X  <  —d 

— V'o,  —d  <  X  <  d 

m4  X  X  +  64,  d  <  X  <  do 


0,  X  >  Go  J 


(3.17) 


where  —do  is  the  boundary  between  Region  I  and  Region  II  in  Figure  3.7,  m2  and  62 
are  the  slope  and  intercept  of  the  trapezoidal  side  in  Region  II,  —a  is  the  boundary 
between  Region  II  and  Region  III,  Vo  is  the  depth  of  the  potential  well  in  Region 
III,  1714  and  64  are  the  slope  and  intercept  of  the  trapezoidal  side  in  Region  IV,  and 
do  is  the  boundary  between  Region  IV  and  Region  V.  Changes  in  the  spatial  grid 
resolution  where  the  potential  is  given  by  Equation  3.17  has  the  effect  of  adding 
more  points  to  the  sloping  sides  of  the  trapezoid  and  not  changing  the  slope  of  the 
trapezoid  to  approach  the  analytical  square  well.  This  ensures  that  the  order  of 
convergence  for  Ax  is  due  solely  to  changes  in  the  spatial  grid  resolution  and  not 
to  the  numerical  square  well  (the  trapezoidal  well)  approaching  the  analytic  square 
well. 

*’See  the  discussion  near  the  end  of  Section  3.1.1. 
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Figure  3.7  One  dimensional  trapezoid  potential  well  used  in  determining  the  order 
of  convergence  for  changes  in  Ax. 

The  efficiency  of  the  numeric  solution,  also  known  as  the  order  of  convergence, 
is  derived  in  the  following  manner.  As  the  resolution  of  the  spatial  grid  is  increased 
by  a  factor  S  =  Axi/Axj+i,  Ax,  >  Ax,^.i  where  Ax  is  the  spatial  grid  resolution, 
the  error  in  the  approximate  solution  should  be  reduced  according  to 


—  =  S’^,  (3.18) 

Ci+l 


where  n  is  the  order  of  convergence.  The  relative  errors  e,  and  e,+i  are  given  by 


— 


(3.19) 


where  Soc  is  the  benchmark  solution.  On  a  log-log  plot  of  the  error  versus  the  grid 
resolution  refinement,  the  order  of  convergence  is  the  slope  of  the  resulting  curve. 
Methods  with  a  large  order  of  convergence  n  show  a  greater  reduction  in  error  for  a 
given  grid  refinement  as  opposed  to  those  methods  with  a  small  order  of  convergence. 
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Figure  3.8  In  a  one  dimensional  trapezoidal  well,  relative  error  vs.  changes  in 
spatial  grid  resolution  for  the  probability  of  transmission  in  a  one  di¬ 
mensional  square  well  potential. 

Figure  3.8  illustrates  the  convergence  behavior  of  the  channel  packet  method 
for  spatial  grid  resolution  refinement.  The  four  runs,  where  for  each  run  the  shape 
of  the  trapezoid  was  changed  in  order  to  reflect  a  trapezoid  approaching  a  square 
well,  exhibit  similar  behavior  and  all  four  have  an  order  of  convergence  of  3.  The 
dramatic  change  and  behavior  in  the  relative  error  for  a  spatial  grid  resolution  of 
greater  than  0.06  au  reflects  the  rules  laid  down  in  Section  3.1.1  for  choosing  the 
grid.  Using  Ax  >  0.06  au  results  in  a  maximum  momentum  grid  value  that  is  too 
small  to  support  the  wave  function  on  the  momentum  grid  which  leads  to  spurious 
behavior  and  incorrect  results  for  the  probability  of  reaction. 

Likewise,  Figure  3.9  illustrates  the  convergence  behavior  of  the  channel  packet 
method  for  temporal  grid  resolution  refinement.  Since  changes  in  At  do  not  affect  the 
shape  of  the  numeric  trapezoidal  well,  only  one  potential  was  needed.  For  changes 
in  the  temporal  grid  resolution,  the  channel  packet  method  exhibits  an  order  of 
convergence  of  4.  .4gain,  the  dramatic  change  in  relative  error  and  the  oscillations 
for  At  >  0.001  au  reflects  a  violation  of  the  rules  in  Section  3.1.1. 
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Figure  3.9  In  a  one  dimensional  square  well,  relative  error  vs.  changes  in  temporal 
grid  resolution  for  the  probability  of  transmission  in  a  one  dimensional 
square  well  potential. 

3.2  Gaussian  well  with  Gaussian  barriers 

A  Gaussian  well  with  Gaussian  barriers  provides  a  smooth  potential  to 
investigate  the  influence  of  barriers  and  wells  on  one  dimensional  scattering.  Since 
the  Gaussian  well  lacks  the  discontinuity  of  the  one  dimensional  square  well,  the 
spatial  grid  requires  fewer  points  in  order  to  effectively  model  the  potential  and  the 
computation  runs  in  less  time.  This  allows  investigation  of  many  different  aspects 
of  one  dimensional  scattering  in  a  relatively  short  time. 

3.2.1  The  potential.  The  potential  is  a  Gaussian  well  with  symmetric 

Gaussian  barriers, 

V  (x)  -  .4exp|-  (a:  -  a:i,o)^}  -  0.09 exp  +  Aexp{-(x  -  xa.o)^}  ,  (3.20) 

where  A  is  the  height  of  the  barriers  and  xi,o  and  X2,o  are  the  left  and  right  offsets 
for  the  barriers  with  the  offsets  chosen  so  that  the  barriers  are  zero  before  the  well 
becomes  non-zero.  Figure  3.10  illustrates  the  potential  for  0.03  au  barriers.  The 
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Table  3.2  Initial  parameters  for  scattering  from  a  one  dimensional  Gaussian  barrier- 
well-barrier  potential  for  a  barrier  height  of  0  au. 


Parameter 

Value  (au) 

Note 

mass 

1224.0 

Reduced  mass  of  H  +  H2  system 

Ax 

X  grid  resolution 

nstep 

number  of  time  steps  for  reactant  propagation 

npstep 

1500 

number  of  time  steps  for  product  propagation 

At 

1.0 

time  step 

ioff 

0.0 

initial  x  offset 

1.5 

initial  momentum  offset  in  units  of 

0.55 

initial  wave  packet  spread  in  coordinate  representation 

initial  parameters  for  the  calculation  are  shown  in  Table  3.2.  The  parameter  imom 
is  changed  with  changes  in  barrier  height  so  that  the  initial  wave  packet  has  a  peak 
energy  near  the  energy  of  the  barrier  height.  Both  the  reactant  and  product  wave 
packets  are  propagated  out  to  the  asymptotic  limit  under  Ho  to  t  =  ±nstep  x  At 
and  back  to  t  =  0  under  the  full  scattering  Hamiltonian  H.  Figure  3.11  illustrates 
the  resulting  reactant  Mpller  state  for  a  barrier  height  of  0.03  au.  The  correlation 
function  is  then  calculated  for  negative  times  from  t  =  0tot  =  — lx  nstep  x  At. 
Since  this  action  reverses  the  propagation  from  the  asymptotic  limit  to  the  interaction 
region,  the  evolving  Mpller  state  quickly  exits  the  interaction  region.  However,  for 
t  >  0  the  evolving  Mpller  state  becomes  trapped  in  the  interaction  region  and  the 
propagation  time  necessary  for  the  correlation  function  to  reach  zero  can  become 
very  long.  The  parameter  npstep  sets  the  limit  for  the  propagation  forward  in 
time  to  t  =  npstep  x  At.  For  a  barrier  height  of  0.05  au,  the  correlation  function 
reached  zero  after  approximately  200, 000  steps  compared  to  the  correlation  function 
reaching  zero  after  1500  time  steps  for  a  barrier  height  of  zero.  Again,  with  non-zero 
barrier  heights  and  without  absorbing  boundary  conditions,  the  grid  would  have  to 
be  quite  large  in  order  to  support  the  evolving  Mpller  state  before  completely  exiting 
the  interaction  region.  These  calculations  were  performed  on  a  spatial  grid  of  2048 
points.  Without  absorbing  boundary  conditions,  the  spatial  grid  would  have  to  be  21 
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times  larger  in  order  to  accommodate  the  evolving  wave  packet  until  the  quasi-bound 
state  completely  exited  the  interaction  region  in  the  case  where  the  barrier  height 
was  0.05  au.  However,  the  restriction  that  the  number  of  points  in  the  spatial  grid 
be  a  pow'er  of  2  forces  the  grid  to  be  32  times  larger  instead  of  21.  An  estimate  can 
be  made  of  the  cost  savings  resulting  from  the  smaller  grid.  Since  the  majority  of  the 
computational  effort  is  in  the  FFTs,  we  can  estimate  the  reducting  in  computational 
time  due  to  the  smaller  grid  from  how  an  FFT  scales  computationally.  For  N  a  power 
of  2,  the  computational  cost  of  an  FFT  scales  as  A^log2  N.  Using  a  small  grid  of  Ap¬ 
points  instead  of  a  larger  grid  of  Ni  points  leads  to  a  savings  in  the  computational 
effort  in  the  FFT  of 


Ni  log2  Ni 
Nj  log2  Nj  ’ 


(3.21) 


where  ^  is  the  estimated  savings.  For  the  case  where  the  barrier  height  was  0.05  au, 
Nj  was  2048  and  A;  is  estimated  to  be  65536.  This  leads  to  an  estimated  savings  of 
^  =  46. 
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Figure  3.11  Reactant  M0ller  state  for  a  Gaussian  barrier- well-barrier  potential. 

The  dashed  line  is  the  potential  and  the  solid  line  is  the  absolute  value 
of  the  Mpller  state,  {x|>P+((  =  0)). 


3.2.2  S-matrix  elements.  Figure  3.12  illustrates  the  probability  of 

transmission  for  the  Gaussian  barrier- well-barrier  potential.  The  Gaussian  well  with¬ 
out  barriers  exhibits  100%  transmission  for  the  range  of  energies  examined  using  the 
channel  packet  method.  The  high  transmission  rate  is  not  unexpected  due  to  the 
smooth  nature  of  the  potential.  The  behavior  at  very  low  energies  is  not  resolvable 
in  the  channel  packet  method  due  to  the  requirement  that  the  initial  wave  packet 
contain  no  negative  momenta.  With  the  addition  of  barriers  at  x  =  ±6  au,  the 
probability  of  transmission  exhibits  significant  amplitude  below  the  barrier  height, 
and  some  reflection  at  energies  above  the  barrier  height.  Transmission  coefficients 
computed  for  all  three  barrier  heights  have  the  same  general  form.  However,  the 
features  in  each  curve  are  offset  in  energy  by  the  barrier  height.  The  influence  of  the 
barriers  is  most  significant  in  the  quasi-bound  states  that  they  trap  in  the  interaction 
region.  A  quasi-bound  state  can  be  pictured  as  the  part  of  the  evolving  Mpller  state 
that  remains  behind  in  the  interaction  region  after  the  major  part  of  the  evolving 
Mpller  state  has  exited  the  interaction  region.  The  remaining  wave  packet  is  not 
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Table  3.3  Initial  parameters  used  for  the  absorbing  boundary  conditions. 


Parameter 

Value 

B 

5.6 

A 

1.0  X  10"' 

Xo 

±7.3 

truly  bound;  over  time  the  quasi-bound  state  will  exit  the  interaction  region.  The 
overlap  between  the  quasi-bound  state  and  the  non-evolving  reactant  Mpller  state 
prevents  the  correlation  function  from  rapidly  approaching  zero.  Without  absorbing 
boundary  conditions,  the  grids  would  have  to  be  quite  large  in  order  to  support  the 
evolving  Moller  state  in  the  presence  of  quasi-bound  states.  For  the  case  where  the 
barrier  height  is  0.03  au,  the  spatial  grid  would  have  to  be  8  times  larger. 

3.2.3  Influence  of  absorbing  boundary  conditions.  The  Gaussian 

barrier-well-barrier  potential,  with  a  barrier  height  of  0.03  au,  was  used  to  inves¬ 
tigate  the  influence  of  absorbing  boundary  conditions  on  the  S-matrix  elements.  In 
particular,  how  do  poorly  chosen  absorbing  boundary  conditions  influence  the  final 
result?  Three  sets  of  poorly  chosen  absorbing  boundary  conditions  were  investigated: 
too  steep,  too  shallow  and  large  initial  step  height.  The  numeric  implementation  of 
absorbing  boundary  conditions,  given  in  Equation  2.42,  leads  to  a  initial  step  where 
the  absorbing  boundary  conditions  are  first  non-zero.  The  Gaussian  potential  was 
chosen  for  two  reasons:  the  overall  computation  time  is  short  allowing  for  a  wide  va¬ 
riety  of  poor  absorbing  boundary  conditions  to  be  investigated,  and  the  probability 
of  transmission  has  many  features  for  poorly  chosen  absorbing  boundary  conditions 
to  influence.  The  initial  parameters  used  for  the  absorbing  boundary  conditions  in 
the  probabilities  of  transmission  illustrated  in  Figure  3.12  are  given  in  Table  3.3. 

Absorbing  boundary  conditions  which  are  too  steep  will  reflect  more  strongly 
than  properly  chosen  absorbing  boundary  conditions.  Table  3.4  lists  the  values  of 
the  parameter  B,  which  governs  the  steepness  in  Equation  2.8,  for  the  three  sets  of 
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Table  3.4  Values  for  the  parameter  B  used  in  determining  when  absorbing  bound¬ 
ary  conditions  are  too  steep. 


Run 

B 

0 

5.6 

1 

2.0 

2 

1.0 

3 

0.1 

4 

0.01 

S-matrix  elements  shown  in  Figure  3.13.  Run  0  is  the  value  of  B  for  the  original  ab¬ 
sorbing  boundary  conditions  that  were  used  in  the  calculations  of  S-matrix  elements 
shown  in  Figure  3.12.  Runs  1  and  2  produced  results  that  were  nearly  indistinguish¬ 
able  from  the  results  of  Run  0.  Runs  3  and  4  produced  the  dramatic  changes  in  the 
S-matrix  elements  shown  in  Figure  3.13.  Absorbing  boundary  conditions  which  are 
too  steep  lead  to  nearly  periodic  oscillations  in  the  S-matrix  elements  that  rapidly 
obscure  the  correct  result.  The  erroneous  probability  of  transmission  elements  range 
from  slightly  more  than  1  to  values  much  more  than  1. 

Absorbing  boundary  conditions  which  are  too  shallow  will  not  fully  attenuate 
the  evolving  Mpller  state.  For  example,  a  wave  packet  headed  out  the  product  chan¬ 
nel,  moving  to  the  right,  will  reappear  in  the  asymptotic  limit  of  reactant  channel,  on 
the  left,  still  moving  to  the  right.  This  is  due  to  the  imposition  of  periodic  boundary 
conditions  by  the  use  of  FFTs.  The  FFT  treats  the  end  points  of  the  grid  as  the 
same  point  thus  the  evolving  wave  packet  wraps  around.  Table  3.5  lists  the  values 
of  the  parameter  5,  which  governs  the  steepness  in  Equation  2.8,  for  the  three  sets 
of  S-matrix  elements  shown  in  figure  3.14.  Again,  Run  0  is  the  original  absorbing 
boundary  condition  computation.  Runs  5  and  6  produced  results  nearly  indistin¬ 
guishable  from  Run  0.  Run  7  leads  to  a  slight,  nearly  periodic  oscillation  overlying 
the  answer  in  Run  0.  Run  8  leads  to  clear  oscillations  which  are  beginning  to  obscure 
the  correct  answer.  While  absorbing  boundary  conditions  which  are  too  shallow  do 
not  fully  attenuate  the  evolving  wave  packet,  they  can  attenuate  a  significant  por- 
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Table  3.5  Values  for  the  parameter  B  used  in  determining  when  absorbing  bound¬ 
ary  conditions  are  too  shallow. 


Run 

B 

0 

5.6 

5 

10.0 

6 

20.0 

7 

40.0 

8 

60.0 

Table  3.6  Values  for  the  step  height  parameter  A  used  in  determining  when  the 
absorbing  boundary  condition  step  height  is  too  high. 


Run 

A 

0 

1  X  10'^ 

9 

1  X  10“^ 

10 

1  X  10"^ 

tion  of  the  evolving  wave  packet  introducing  less  error  into  the  S-matrix  elements 
than  absorbing  boundary  conditions  which  are  too  steep.  Figure  3.14  illustrates  the 
results  when  the  absorbing  boundary  conditions  are  too  shallow. 

The  final  set  of  poorly  chosen  absorbing  boundary  conditions  concerns  the 
initial  step  height  governed  by  the  parameter  A  in  Equation  2.8.  Initially,  the  step 
height  was  fixed  at  the  very  low  value  of  1x10“^  in  an  effort  to  avoid  any  reflection 
from  a  sudden  onset  of  absorbing  boundary  conditions.  Table  3.6  lists  the  values 
of  the  parameter  A  for  the  three  sets  of  S-matrix  elements  shown  in  Figure  3.15. 
The  step  height  was  adjusted  as  well  as  the  steepness  B  so  that  the  absorbing 
boundary  conditions  had  the  same  peak  amplitude  as  in  Run  0.  Again,  Run  0 
is  the  original  result.  As  seen  in  Figure  3.15,  the  step  height  has  only  a  slight 
influence  on  the  final  values  of  the  S-matrix  elements.  From  these  computations,  a 
valid  first  choice  of  absorbing  boundary  conditions  has  a  parameter  A  of  1  x  10“^ 
and  a  parameter  B  such  that  the  absorbing  boundary  conditions  are  as  shallow  as 
possible  while  still  reaching  a  peak  amplitude  greater  than  7  au.  A  valid  initial  set  of 
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absorbing  boundary  condition  parameters  based  on  this  analysis  is  A  =  1  x  10“^  and 
B  =  5.6.  The  value  of  xq  is  dependent  on  the  reaction  under  consideration  can  not 
be  determined  in  advance  since  it  sets  the  point  at  which  the  absorbing  boundary 
conditions  are  non- zero. 

3.3  Summary 

Since  S-matrix  elements  for  the  one  dimensional  square  well  have  an  ana¬ 
lytic  solution,  a  comparison  of  the  analytic  solution  with  the  computational  result, 
illustrated  in  Figure  3.5,  shows  the  combination  of  the  channel  packet  method  with 
absorbing  boundary  conditions  is  a  valid  approach  to  computing  S-matrix  elements. 
The  validity  of  the  method  is  further  demonstrated  through  two  tests  of  convergence. 
The  first,  illustrated  in  Figure  3.6,  shows  that  as  the  spatial  grid  resolution  is  refined, 
the  numeric  solution  approaches  the  analytic  solution.  Second,  the  channel  packet 
method  exhibits  third  order  convergence  for  spatial  grid  resolution  refinement  and 
fourth  order  convergence  for  temporal  grid  resolution  refinement.  These  orders  of 
convergence  are  illustrated  in  Figures  3.8  and  3.9. 

Using  a  Gaussian  well  with  symmetric  Gaussian  barriers,  the  effects  of  barriers 
and  wells  on  scattering  was  explored.  A  well  without  barriers  may  or  may  not  trap 
a  quasi-bound  state  in  the  interaction  region  requiring  long  computation  times  in 
calculating  the  correlation  function.  However,  the  addition  of  barriers  does  lead 
to  quasi-bound  states  with  possibly  long  lifetimes.  Correct  S-matrix  elements  are 
rapidly  obscured  by  errors  introduced  into  the  correlation  function  calculation  by 
poorly  chosen  absorbing  boundary  conditions  as  seen  in  Section  3.2.3.  Absorbing 
boundary  conditions  which  are  too  steep  quickly  lead  to  large  oscillations  in  the 
probability  of  transmission  and  unphysical  values  much  greater  than  one.  Absorbing 
boundary  conditions  which  are  too  shallow  lead  to  oscillations  in  the  probability  of 
transmission  that  rapidly  obscure  the  correct  result.  However,  the  initial  step  height 
introduced  by  the  numeric  implementation  of  absorbing  boundary  conditions  has 
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only  a  small  influence  on  the  probability  of  transmission  and  a  choice  of  1  x  10  ^  for 
the  parameter  A  in  Equation  2.42  is  a  valid  first  choice. 
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IV.  Two  dimensional  H  +  H2 

The  two  dimensional  H  +  H2  system  has  been  extensively  studied  and 
S-matrix  elements  have  been  calculated  using  a  wide  variety  of  methods. 

In  this  study,  computing  S-matrix  elements  for  collinear  H  +  H2  and  comparing  the 
results  with  previous  work  provides  proof  of  concept  for  the  combination  of  absorbing 
boundary  conditions  with  the  channel  packet  method.'*^ 

4.1  Background 

Liu,  Siegbahn,  Truhlar  and  Horowitz  (LSTH)  have  computed  a  highly  ac¬ 
curate  potential  energy  surface  for  collinear  H  +  The  LSTH  potential  energy 

surface,  which  is  accurate  to  within  0.1%  of  the  true  potential,  is  characterized  by 
two  troughs  parallel  to  the  relative  coordinate  in  each  of  the  channels  H  +  H2  and 
H2  +  H,  and  a  saddle  point,  or  barrier,  in  the  interaction  region.  Computations 
involving  the  H  +  H2  system  are  simplified  by  the  barrier  in  the  interaction  region. 
Figure  4.1  is  a  contour  plot  of  the  LSTH  potential  energy  surface.  The  plain  in  the 
upper  right  corner  of  Figure  4.1  represents  the  H+H+H  channel  and  is  energetically 
inaccessible  in  this  study. 

The  application  of  the  channel  packet  method  begins  with  the  definition  of 
two  initial  wave  packets  (X,  K)  and  {^1 T)  at  f  =  0,  where  X  and  Y 
are  the  bond  coordinates  describing  the  three  collinear  hydrogen  atoms. The  first 
wave  packet  used  to  describe  the  incoming  reactants,  H  H2  (v),  in 

arrangement  channel  1  where  the  diatom  is  prepared  in  a  single  vibrational  eigenstate 
labeled  u.  A  Gaussian  wave  packet  is  used  to  describe  the  wave  packet  in  the 
relative  coordinate  between  H  and  H2.  The  second  w'ave  packet  (X,  T)  is 
used  to  represent  outgoing  products  H2  (V)  +  H  in  arrangement  channel  2  where 
the  diatom  is  in  a  single  vibrational  eigenstate  labeled  V.  Again,  a  Gaussian  wave 
packet  is  used  to  describe  the  wave  packet  in  the  relative  coordinate  between  JI2 
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and  H.  The  next  step  is  to  compute  a  pair  of  intermediate  states  obtained  from 
the  initial  states  by  propagating  (A',  V')  backward  in  time  to  t  =  -oo  using  the 

asymptotic  channel  Hamiltonian  and  propagating  ^’)  forward  in  time 

to  t  =  +CXD  using  the  channel  Hamiltonian  Hq.  This  is  illustrated  in  Figure  4.2 
where  infinity  is  reached  numerically  at  t  =  -4000  au.  The  intermediate  reactant 
(product)  states  are  then  propagated  from  +oo(— oo)  using  the  full  Hamiltonian 
backward  (forward)  to  time  ^  =  0  to  obtain  the  reactant  (product)  Mpller  states. 
The  correlation  function  between  the  resulting  reactant  and  product  Mpller  states 
is  then  computed  using  absorbing  boundary  conditions  as  they  continue  to  evolve 
forward  and  backw'ard  in  time  to  ^  =  ioo.  The  application  of  absorbing  boundary 
conditions  wdth  the  channel  packet  method  when  computing  the  correlation  function 
is  facilitated  by  the  fact  that  the  resulting  reactant  and  product  Mpller  states  are 
typically  well  localized  in  the  interaction  region  of  the  potential  at  t  =  0.  The  well 
localized  Moller  states  (AT,  V’)and’I'i’°  (A'^,  T)  are  shown  in  Figures  4.3  and  4.4. 
To  compute  the  correlation  function,  we  choose  one  M0ller  state  to  evolve  while  the 
other  remains  static.  The  evolving  Mpller  state  will  eventually  reach  the  absorbing 
boundary  conditions  and  be  exponentially  attenuated.  However,  since  the  absorbing 
boundary  conditions  are  placed  where  the  static  Mpller  state  has  zero  amplitude, 
there  is  no  overlap  between  the  static  Mpller  state  and  the  evolving  Mpller  state  in 
the  region  where  the  evolving  state  is  being  attenuated.  The  correlation  function 
is  therefore  unaffected  by  the  absorbing  boundary  conditions.  Figure  4.5  illustrates 
the  placement  of  absorbing  boundary  conditions  with  respect  to  the  Mpller  states. 
S-matrix  elements  are  then  obtained  from  the  Fourier  transform  of  the  correlation 
function. 

Recently,  Jackel  and  Meyer'*®  employed  absorbing  boundary  conditions  in  a 
version  of  the  channel  packet  method  adapted  to  work  with  a  time-dependent-self- 
consistent-field  approach  to  compute  S-matrix  elements  for  the  collinear  H  +  H2 
reaction.  Dai  and  Zhang^®  have  also  used  absorbing  boundary  conditions  with  the 
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Figure  4.5  Contour  plot  of  the  placement  of  the  absorbing  boundary  conditions 
with  respect  to  the  Mpller  states.  The  two  closed  contours  are  the  0.1 
au  contours  of  the  Mpller  states.  The  two  parallel  contours  are  the  0.1 
au  contours  of  the  LSTH  potential  energy  surface.  The  saddle  point  is 
marked  with  a  dot.  The  shaded  region  represents  the  area  in  which  the 
absorbing  boundary  conditions  are  non-zero. 


channel  packet  method  in  their  analysis  of  the  three  dimensional  H  +  O2  reaction. 
In  their  analysis,  Dai  and  Zhang  selected  a  product  Mpller  state  located  in  the 
asymptotic  limit.  As  a  result,  the  absorbing  boundary  conditions  are  placed  far 
from  the  interaction  region  of  the  potential.  The  idea  behind  this  research  is  to  use 
the  smallest  possible  grid  in  order  to  achieve  the  most  efficient  use  of  computational 
resources.  Any  grid  must  contain  the  interaction  region  where  the  physics  occurs. 
A  small  grid  must  somehow  eliminate  the  asymptotic  regions  while  simultaneously 
supporting  the  evolving  wave  packet  until  the  interaction  is  complete.  A  product 
Mpller  state  that  lies  in  the  asymptotic  region  far  from  the  interaction  region  leads 
to  the  use  of  a  large  grid  which  defeats  the  purpose.  Our  approach  differs  in  that  we 
localize  both  the  reactant  and  product  Mpller  states  in  the  interaction  region.  Thus 
there  is  no  need  to  support  the  far  regions  of  the  asymptotic  limit,  w'here  Dai  and 
Zhang’s  product  Mpller  state  is  placed,  leading  to  the  use  of  a  smaller,  more  efficient 
grid. 

4-2  Computational  Procedure 

The  initial  wave  packets  are  constructed  from  a  linear  combination  of 
eigenstates  of  the  relative  Hamiltonian  Hreh  given  by  Equation  2.21  and  a  single 
eigenstate  of  the  internal  Hamiltonian  Hint-  The  parameters  deltaxl,  deltax2,  and 
deltat  are  all  chosen  according  to  Kosloff’s  paper®®  and  discussed  in  more  detail  in 
Section  3.1.1.  The  choice  of  xminl  and  xmin2  is  based  on  the  behavior  of  the  LSTH 
potential  near  zero.  Near  zero,  the  LSTH  potential  flattens  out  at  a  value  of  5  au 
and  xminl  and  xmin2  are  chosen  so  that  the  potential  at  {xminl, xmin2)  is  just  at 
5  au.  Since  the  time  necessary  to  fully  compute  the  correlation  function  is  greatly 
affected  by  the  grid  size,  compact  Mpller  states  are  necessary  in  order  to  use  small 
grids.  In  the  one  dimensional  case,  the  final  form  of  the  Mqller  state  was  heavily 
influenced  by  the  parameters  ioff,  imom  and  sprd.  Likewise,  the  final  form  of  the 
reactant  Mpller  state  in  channel  1  is  heavily  influence  by  ioffl,  imoml  and  sprdl. 
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Following  the  results  of  the  one  dimensional  square  well  case,  the  initial  choice  for 
ioffl  was  made  to  center  the  initial  wave  packet  over  the  interaction  region  with  a 
wave  packet  spread,  sprdl,  such  that  the  initial  wave  packet  did  not  extend  past 
the  edge  of  the  interaction  region.  The  choice  of  the  momentum  offset,  imoml,  was 
based  on  the  requirement  that  the  initial  wave  packet  have  no  negative  momentum 
in  the  relative  coordinate.  The  reactant  Moller  state  was  computed  based  on  these 
choices.  However,  this  turned  out  to  not  be  the  optimal  choice  for  two  reasons. 

In  order  to  accommodate  the  spread  of  the  wave  packet  on  the  spatial  grid  as  the 
packet  is  propagated  back  in  time,  a  long  narrow  grid  is  the  most  efficient  one  to  use. 
The  spread  along  with  the  long  narrow  grid  is  illustrated  in  Figure  4.2.  In  centering 
the  initial  wave  packet  in  the  interaction  region,  the  resulting  Mpller  state  at  t  =  0 
has  bifurcated  into  both  channels  and  the  initial,  long,  narrow  grid  would  have  to 
be  widened  to  accommodate  the  bifurcation.  Instead,  the  parameter  ioffl  is  set  so 
that  the  initial  wave  packet  is  slightly  offset  into  the  reactant  channel.  The  resulting 
Moller  state  is  well  localized  in  the  interaction  region  as  illustrated  in  Figure  4.3. 
The  parameters  for  the  reactant  Mpller  state  calculation  are  shown  m  Table  4.1.  The 
H+H2  potential  energy  surface  has  a  barrier  and  no  well  so  the  total  energy  contains 
only  the  maximum  relative  kinetic  energy  and  internal  vibrational  eigenstate  energy. 

To  obtain  eigenstates  of  the  internal  Hamiltonian,  Hint  is  represented  in  a  ba¬ 
sis  of  Morse  oscillator  eigenstates  and  that  matrix  representation  is  diagonalized. 
The  transformation  matrix  that  diagonalizes  Hint  contains  the  linear  combination 
of  Morse  oscillator  eigenstates  required  to  compute  the  eigenstates  of  Hint-  In  the 
collinear  H  +  H2  reaction,  the  Mpller  state  calculations  are  most  efficiently  com¬ 
puted  on  long,  rectangular  grids  as  illustrated  in  Figure  4.2.  The  initial  wave  packet 
is  shown  at  t  =  0  with  the  intermediate  state  at  t  =  -4000  au.  The  intermediate 
state  is  obtained  by  analytically  propagating  the  initial  state  from  t  =  0tot  =  -4000 
au  using  Jacobi  coordinates.  The  propagation  spreads  the  initial  wave  packet  m  the 
relative  coordinate  r  according  to  Equation  2.34  and  develops  a  phase  factor  in  the 
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Initial  Mciller  state  calculation  parameters  for  the  H  +  H2  reactant 
channel. 


Parameter 

Value  (au) 

Notes 

xminl 

1x10-' 

minimum  X  grid  boundary 

xmin2 

1x10"' 

minimum  Y  grid  boundary 

ma 

1837.1526025 

mass  of  hydrogen  atom  a 

mb 

1837.1526025 

mass  of  hydrogen  atom  b 

me 

1837.1526025 

mass  of  hydrogen  atom  c 

deltaxl 

0.2 

X  grid  resolution 

deltax2 

0.2 

Y  grid  resolution 

nstep 

4000 

number  of  time  steps 

sftent 

4000 

number  of  time  steps  to  stop  the  propagation  at 

deltat 

1.0 

time  step 

ioffl 

1.40105 

initial  X  offset 

iofF2 

0.0 

initial  Y  offset 

imoml 

-1.3 

initial  Px  offset 

imom2 

0.0 

initial  Py  offset 

sprdl 

0.4 

initial  relative  wave  packet  spread  in  the  X  direction 

sprd2 

0.0 

initial  relative  wave  packet  spread  in  the  Y  direction 

temporal  offset  for  channel  1  initial  wave  packet 

temporal  offset  for  channel  2  initial  wave  packet 

evch 

1 

asymptotic  eigenstate  number 
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internal  coordinate  R  according  to  Equation  2.33.  The  intermediate  wave  packet 
is  then  transformed  to  bond  coordinates  and  the  final  Moller  state  is  computed  by 
propagating  back  to  t  =  0  under  the  full  scattering  Hamiltonian  H.  This  propagation 
is  done  numerically  using  the  split-operator  method  discussed  in  Section  2.6,  using 
a  modified  version  of  the  code  first  developed  by  Weeks  and  Tannor.^^  Figures  4.3 
and  4.4  show  the  Mpller  states  (A',  V'|'5+°)  and  (A',  respectively.  Note  that 

the  Mailer  states  are  well  localized  in  the  interaction  region.  This  permits  the  ab¬ 
sorbing  boundary  conditions  to  be  placed  close  to  the  interaction  region  and  still 
have  the  spatial  extent  to  fully  attenuate  the  evolving  Mailer  state  during  the  corre¬ 
lation  function  computation.  The  placement  of  the  absorbing  boundary  conditions 
with  the  LSTH  potential  energy  surface  and  the  Mailer  states  is  shown  in  Figure  4.5. 
The  absorbing  boundary  conditions  in  two  dimensions  take  a  form  similar  to  those 
in  one  dimension  given  by  Equation  2.42.  For  the  collinear  /f  -f  i/2  reaction,  the 
absorbing  boundary  conditions  are  given  by 


Va(x,y) 


^  ±h  (Xo)  Ai  exp  { } ,  a:  >  y  ' 
^  ±h  (yo)  Ai  exp  {  ,  X  <  y 


(4.1) 


where  h  (if)  is  the  heavyside  step  function  given  in  Equation  2.43,  Xq  and  Vo  are  where 
the  absorbing  boundary  conditions  begin  in  channel  1  and  channel  2,  respectively, 
and  A  and  B  govern  the  shape  of  the  absorbing  boundary  conditions.  For  the 
collinear  H  -\-  H2  reaction,  the  parameters  A  =  1  x  10“^  and  B  =  2  were  chosen  so 
that  the  evolving  Mpller  state  was  completely  attenuated  before  reflecting  from  the 
grid  edges  with  a  minimum  of  reflection.  The  offsets  Xq  and  Yq  were  chosen  to  be 
6.2  so  that  the  absorbing  boundary  conditions  begin  just  beyond  where  the  Mpller 
states  reached  numeric  zero. 

In  Weeks  and  Tannor’s  original  calculation,^’  the  grid  size  was  256  points  by 
256  points  in  order  for  the  correlation  function  to  reach  zero  before  the  evolving 
Mpller  state  reached  the  edge  of  the  grid  and  reflected.  The  reactant  and  product 
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Moller  states  are  well  localized  inside  a  grid  of  32  by  32  points.  An  initial  grid  of  64 
by  64  points  was  used  to  apply  the  combination  of  absorbing  boundary  conditions 
and  the  channel  packet  method  to  the  computation  of  the  correlation  function.  The 
choice  of  64  x  64  was  partly  arbitrary  and  partly  due  to  the  FFT  being  most  efficient 
for  :V  a  pow'er  of  2. 

4-3  Results 

In  order  to  compute  S-matrix  elements  for  the  i/  =  0  — )•  i/'  =  0, 1  reactions, 
the  Moller  state  is  computed  using  the  =  0  internal  vibrational  eigenstate 

and  two  Mpller  states  are  computed  using  the  t/'  =  0,1  internal  vibrational 

eigenstates.^^  At  t  =  0,  the  reactant  and  product  Mpller  states  are  localized  in  the 
interaction  region  of  the  potential  and  the  absorbing  boundary  conditions  are  located 
as  shown  in  Figure  4.5.  In  this  computation,  the  reactant  Mpller  state  is  chosen  to 
evolve  in  time  to  both  t  =  —t  and  t  =  +r.  As  the  reactant  Mpller  state  evolves 
backward  in  time,  it  retraces  its  path  along  channel  1.  However,  as  it  reaches  the 
region  containing  the  absorbing  boundary  conditions,  it  is  steadily  attenuated  before 
it  reaches  the  edge  of  the  smaller  grid.  While  there  is  some  reflection,  it  is  minimized 
by  the  choice  of  the  absorbing  boundary  conditions.  Likewise,  as  the  reactant  Mpller 
state  evolves  forward  in  time,  it  bifurcates  into  both  channels.  .Again,  as  it  reaches 
the  regions  containing  the  absorbing  boundary  conditions,  it  is  attenuated  with  a 
minimum  of  reflection  and  does  not  reach  the  grid  edges  as  illustrated  in  Figure  4.6. 
As  the  reactant  Mpller  state  evolves  from  — r  to  r,  its  correlation  Cy-y  {t)  with  the 
product  Mpller  state  is  calculated.  Correlation  functions  for  the  i/  =  0  to  u'  =  0 
and  1/  =  0  to  i/'  =  1  reactions  are  shown  in  Figures  4.7  and  4.8.  The  correlation 
function  C-y'-y  (t)  is  then  used  to  compute  S-matrix  elements,  5^/^  ,  according  to 
Equation  2.30.  Figures  4.9  and  4.10  show  the  probability  of  transmission  for  two 
reactions  H  +  H2  {i'  =  0)  H2{v''  —  0,1)  +  i/.and  are  in  agreement  with  previous 

calculations. In  order  to  quantify  the  error  in  the  new  method  compared 
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Figure  4.6  Contour  plot  of  the  time  evolution  of  the  Mpller  state  at  t  =  +1000 

au.  The  closed  contours  are  the  0.1  au  contours  of  |(X, yi^+°)|.  The 
thick  closed  contour  is  the  0.1  au  contour  of  the  non-evolving  Mpller 
state  |(X,  y  |^?.’°)|.  The  two  parallel  contours  are  the  0.1  au  contours  of 
the  LSTH  potential  energy  surface.  The  shaded  region  represents  the 
area  in  which  the  absorbing  boundary  conditions  are  non-zero. 


Figure  4.7  Absolute  value  of  the  correlation  function  calculated  using  the 
channel  packet  method  with  absorbing  boundary  conditions  for 
H+H2{v  =  0)-^H2{v’=0)+H. 


Figure  4.8  Absolute  value  of  the  correlation  function  calculated  using  the 
channel  packet  method  with  absorbing  boundary  conditions  for 
H  +H2(v  =  0)  ^  H2(v’=  1)+H  . 
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Figure  4.9  Probability  of  reaction  calculated  using  the  channel  packet 
method  with  absorbing  boundary  conditions  for 
if  +H2(v  =  0)  -»  if2(v  =  0)  +  ii. 


Energy  (eV) 


Figure  4.10  Probability  of  reaction  calculated  using  the  channel  packet 
method  with  absorbing  boundary  conditions  for 

H+H2{v  =  0)^H2{v^=1)  +  H. 
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to  previous  work,  the  work  of  Weeks  and  Tannor^^  was  used  as  a  benchmark  and 
the  L\  norm  was  computed  using 


L, 


1  \Pweeks  ^a6c| 
N  ^  Pweeks 


(4.2) 


where  N  is  the  number  of  probabilities  of  reaction  computed,  pweeks  is  the  work  of 
Weeks  and  Tannor^®  and  pabc  is  the  computation  performed  using  absorbing  bound¬ 
ary  conditions.  The  L\  norm  error  was  6.3  x  10~^  which  is  at  the  third  digit  conver¬ 
gence  test  used  in  Chapter  III  and  the  calculation  is  considered  to  have  converged 
to  the  correct  solution. 

Using  absorbing  boundary  conditions,  the  correlation  functions  used  to  com¬ 
pute  S-matrix  elements  shown  in  Figures  4.9  and  4.10  are  computed  using  a  64  x  64 
grid  compared  with  a  256  x  256  grid  required  to  compute  the  correlation  function 
without  absorbing  boundary  conditions.^'’  The  calculation  of  without  absorbing 
boundary  conditions  took  ~  31,000  seconds  of  CPU  time  on  a  MIPS  R3000  Sili¬ 
con  Graphics  workstation  and  on  the  same  computer  the  calculation  with  absorbing 
boundary  conditions  took  ~  2500  seconds  of  CPU  time.  Similarly,  on  a  MIPS  RIOOOO 
Silicon  Graphics  workstation,  the  calculation  without  absorbing  boundary  conditions 
took  ~  1450  seconds  of  CPU  time  and  the  calculation  with  absorbing  boundary  con¬ 
ditions  took  ~  80  seconds.  Thus,  an  order  of  magnitude  improvement  in  the  calcu¬ 
lation  time  is  obtained  when  using  absorbing  boundary  conditions  with  the  channel 
packet  method  to  compute  the  correlation  function  for  the  H  +  H2  reaction.  Similar 
to  Chapter  III,  an  estimate  of  the  savings  can  be  computed  from  the  computational 
scaling  of  the  FFTs  used  in  the  calculation.  For  a  grid  of  Ni  x  N2  points,  the  FFT 
scales  as  (iVi  x  N2)  log2  {Ni  x  N2).  In  the  special  case  where  Ni  =  N2,  this  reduces 
to  2  X  N'^  log2  N.  Reducing  the  grid  from  256  x  256  to  64  x  64  leads  to  an  estimated 
cost  savings  of  21.  The  actual  savings  was  18  with  the  difference  due  to  the  fact  that 
not  all  the  computational  effort  is  in  the  FFTs. 
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4-4  Convergence  testing 

Similar  to  the  one  dimensional  case,  the  channel  packet  method  in  two 
dimensions  is  tested  against  two  criteria  for  convergence.  The  first  test  is  convergence 
to  the  correct  answer  and  the  second  is  the  efficiency,  or  order,  of  convergence. 
Convergence  to  the  correct  answer  is  is  necessary  as  a  method  which  converges, 
but  not  to  the  correct  answer,  is  of  little  use.  Figure  4.11  illustrates  the  channel 
packet  method  with  absorbing  boundary  conditions  converging  towards  the  correct 
answer  for  the  reaction  FT  +  i/2  =  0)  ^  ^2  =  0)  +  //.  The  data  set  labels  in 

Figure  4.11  do  not  reflect  the  values  of  any  one  parameter.  Between  each  of  the  data 
sets,  Axi,  Aa:2,  N  and  At  may  have  changed.  The  line  labeled  “10”  in  Figure  4.11 
illustrates  the  result  when  the  guidance  for  grid  selection  in  Section  3.1.1  is  violated. 
In  this  case,  the  temporal  grid  was  too  coarse  to  properly  support  the  wave  packet. 
As  in  Section  3.1.4,  the  Li  norm  defined  in  Equation  3.16  is  used  to  judge  the 
convergence  towards  the  correct  solution.  The  data  sets  labeled  “1”  and  “2”  have 
in  fact  converged  to  the  correct  solution  such  that  the  Li  norm  has  converged  to  3 
digits. 

The  second  test  of  convergence,  the  efficiency  or  order,  was  outlined  in  Sec¬ 
tion  3.1.4.  In  the  two  dimensional  case,  only  the  parameter  At  can  be  isolated  and 
tested  due  to  the  nature  of  the  H  +  H2  collinear  reaction.  The  symmetric  potential 
energy  surface  and  the  two  reaction  channels  require  equal  spatial  grid  parameters. 
Unequal  spatial  grid  parameters  in  the  Mpller  state  calculation  requires  interpolation 
between  grid  points  for  one  or  both  Mpller  states  in  the  correlation  function  calcula¬ 
tion.  The  order  of  convergence  for  the  two  dimensional  channel  packet  method  is  4 
as  is  illustrated  in  Figure  4.12.  This  is  the  same  as  in  the  one  dimensional  case.  The 
dramatic  rise  in  the  relative  error  in  Figure  4.12  for  values  of  At  >  1.0  au  is  due  to 
temporal  grid  that  is  too  coarse  to  support  the  wave  packet. 
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Probability  of  Reaction 


Figure  4. 11  Convergence  of  the  probability  of  reaction  using  the  channel  packet 
method  with  absorbing  boundary  conditions  towards  the  correct  answer  for 
the  reaction  //  +  f/2  ( v  =  0)  — >  //2 ^  • 


Figure  4.12  Relative  error  vs.  Temporal  grid  resolution  for  the  probability  of  reaction 
for  the  reaction  H  +  H2{v  =  0)-^  O)  +  /f . 
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4-5  Summary 

The  two  dimensional  H  +  H2  system  has  been  extensively  studied  and  S- 
matrix  elements  have  been  calculated  using  a  wide  variety  of  methods*'’^'^’^'"^®’^^  and 
provides  an  excellent  benchmark  for  testing  the  combination  of  the  channel  packet 
method  with  absorbing  boundary  conditions.  The  results  are  in  agreement  with 
previous  calculations^^’ with  the  Li  norm  error  between  the  channel  packet 
method  and  the  channel  packet  method  with  absorbing  boundary  conditions  is  on 
the  order  of  6  x  10“^.  This  is  a  small  error  in  comparison  to  the  dramatic  reduction  in 
computational  effort  the  combination  of  the  channel  packet  method  with  absorbing 
boundary  conditions  produces.  The  grid  for  computing  the  correlation  function  was 
reduced  by  a  factor  of  16,  from  256  x  256  to  64  x  64.  Overall,  the  result  of  the  grid 
reduction  was  an  order  of  magnitude  improvement  in  the  time  necessary  to  compute 
the  correlation  function."*^ 
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V.  Two  dimensional  reactive  quantum  scattering  for  a  model 
system  of  two  coupled  Morse  oscillators 

A  simple  model  system  of  two  coupled  Morse  oscillators  is  proposed  to 
investigate  the  combination  of  the  channel  packet  method  with  absorbing  boundary 
conditions  in  systems  with  a  well  in  the  interaction  region.  In  this  chapter,  three 
mass  configurations  are  investigated:  light-light-light  (LLL),  medium-light-medium 
(MLM)  and  heavy-light-heavy  (HLH).  The  light-heavy-light  (LHL)  mass  configura¬ 
tion  is  treated  in  the  next  chapter.  The  mass  configurations  label  the  masses  of  the 
collinear  atoms  in  the  reaction  from  left  to  right.  For  example,  HLH  labels  atom  A 
as  heavy,  B  as  light  and  C  as  heavy  in  the  reaction  A  -f  BC  (r)  — >■  AB  {u')  -p  C . 

5.1  Two  coupled  Morse  oscillator  potential 

The  equation  for  the  potential  energy  surface  of  two  coupled  Morse  oscil¬ 
lators  is  simply 

V  (X,  F)  =  De  {1  -  exp  [-a{X  -  re)]}'  +  De  {1  -  exp  [-a(F  -  r,)]f  -  D,,  (5.1) 

where  is  the  dissociation  energy,  a  is  related  to  the  anharmonicity  and  Tg  is 
the  diatomic  equilibrium  separation.  The  dissociation  energy  is  subtracted  from  the 
potential  energy  surface  so  that  in  the  asymptotic  limit  at  the  equilibrium  position 
the  potential  energy  surface  is  equal  to  zero.  The  parameters  governing  the  potential 
energy  surface  are  completely  variable  in  this  model.  For  the  reactions  in  this  chapter, 
a  single  potential  energy  surface  was  constructed  having  a  dissociation  energy  of  Dg  = 
0.1  au  =  2.72  eV,  an  equilibrium  separation  of  Tg  =  0.7  au,  and  the  anharmonicity, 
a  =  5,  such  that  the  asymptotic  limit  could  be  reached  while  remaining  on  a  256x256 
grid.  These  choices  are  governed  by  the  fact  that  in  order  to  reach  the  asymptotic 
limit,  the  value  of  the  AB  Morse  oscillator  potential  in  the  A  +  BC  channel  has  to  be 
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numerically  zero.  However,  there  must  still  be  room  on  the  coordinate  grid  in  that 
channel  so  that  the  evolving  wave  packets  completely  exit  the  interaction  region, 
i.e.  reach  the  asymptotic  limit,  without  reflecting  from  the  edge  of  the  grid.  For 
comparison,  a  model  proton  transfer  system®*  has  =  0.16  au  and  the  collinear 
H  +  H2  reaction  has  =  0.17  au.  A  contour  plot  of  the  model  potential  energy 
surface  where  Z)e  =  0.1  au  and  a  =  5  is  shown  in  Figure  5.1.  The  well  is  0.1  au 
below  the  bottom  of  the  asymptotic  limit  of  the  troughs  and  0.2  au  below  the  plain 
representing  the  energetically  inaccessible  A  +  B  +  C  channel. 


5.2  iMass  configurations 


Initially,  three  mass  configurations  were  considered:  LLL,  MLM  and  HLH. 
The  purpose  was  to  investigate  for  the  same  potential  energy  surface  the  effect  of 
the  kinetic  energy  operator  T  on  reactive  quantum  scattering.  In  bond  coordinates, 
T  is  given  by 


f  = 


H  I  H 


PxPy 

rrii, 


(5.2) 


where  Px  and  Py  are  the  conjugate  momentum  operators  for  the  bond  coordinates 
(X,  Y)  and  the  reduced  masses  ^aM  and  Hc^ab  are  given  by 


mj  {mk  +  m/) 
mj  +  mfc  +  m/  ’ 


(5.3) 


Equation  5.2  illustrates  that  the  kinetic  energy  operator  is  parameterized  by  the 
masses  of  the  reacting  particles.  In  this  investigation,  the  masses  were  chosen  to 
model  the  H  +  H2  reaction  for  the  LLL  mass  configuration.  The  MLM  configuration 
was  chosen  to  be  slightly  mass  disparate  compared  to  the  non-mass  disparate  LLL 
configuration  and  the  HLH  configuration  was  chosen  to  have  very  disparate  masses. 
The  mass  configurations  are  given  in  Table  5.1.  The  three  configurations  have  differ¬ 
ent  asymptotic  limits  in  the  number  of  internal  vibrational  eigenstates  and  energies 
due  to  the  mass  in  the  kinetic  energy  term  in  Equation  5.2.  The  eigenstates  and 
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Li 


Contour  plot  for  the  model  two  coupled  Morse  oscillator  potential  energy 
surface.  The  contours  range  from  -0. 1  au  to  0. 1  au  in  steps  of  0.04  au.  The 
dot  in  the  center  of  the  well  represents  the  lowest  point  on  the  potential 
energy  surface. 


Figure  5.2  Reactant  M0ller  state  for  the  two  coupled  Morse  oscillator  LLL  mass 
configuration.  The  two  parallel  contours  are  the  0.1  au  contours  of 
the  two  coupled  Morse  oscillator  potential  energy  surface.  The  closed 
contours  represent  the  absolute  value  of  the  Mpller  state  {X, 


plitude  of  (.V,  has  exited  the  interaction  region  into  the  product  channel.  A 

small  amount  of  reflection  back  into  the  reactant  channel  is  also  visible.  The  correla¬ 
tion  function  calculation,  shown  in  Figure  5.5,  rapidly  approaches  zero  for  negative 
times  but  not  for  positive  times  indicating  the  presence  of  a  quasi-bound  state  in  the 
interaction  region. 

The  probability  of  reaction  shown  in  Figure  5.6  exhibits  an  overall  oscillation 
that  is  questionable  in  nature.  Similar  oscillations  can  be  induced  by  terminating 
the  computation  of  the  correlation  function  prior  to  the  correlation  function  reach¬ 
ing  zero.  However  in  Figure  5.5  it  is  clear  that  the  correlation  function  has  reached 
zero  and  the  oscillation  is  not  due  to  early  termination  of  the  correlation  function 
computation.  In  an  attempt  to  understand  the  oscillations  in  Figure  5.6,  the  cor¬ 
relation  function  was  computed  using  a  smaller  grid  where  the  absorbing  boundary 
conditions  have  been  placed  closer  to  the  interaction  region.  The  correlation  func¬ 
tion  for  this  calculation  is  shown  in  Figure  5.7.  Figure  5.8  illustrates  the  probability 
of  reaction  for  the  LLL  mass  configuration  computed  from  the  correlation  function 
shown  in  Figure  5.7.  While  the  underlying  shape  of  the  curve  has  not  changed, 
including  the  drop  in  the  probability  of  reaction  at  0.75  eV,  the  oscillations  have 
changed  in  both  frequency  and  amplitude.  In  order  to  determine  if  the  feature  at 
0.75  eV  is  real  or  an  artifact  of  the  absorbing  boundary  condition  reflection  error, 
the  correlation  function,  shown  in  Figure  5.9,  for  the  LLL  mass  configuration  was 
computed  without  using  absorbing  boundary  conditions.  The  feature  at  0.75  eV  is 

Table  5.2  The  three  lowest  internal  vibrational  eigenstate  energies  for  the  two  cou¬ 
pled  Morse  oscillator  potential  energy  surface  for  the  three  mass  config¬ 
urations.  Energies  are  in  atomic  units. 


M 

LLL 

MLM 

HLH 

El 

0.03223 

0.02967 

0.02551 

1 

0.07794 

0.07339 

0.06528 

2 

0.09865 

0.09627 

0.09006 
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Table  5.3  Initial  M0ller  state  calculation  parameters  for  the  LLL  mass 
configuration. 


Parameter 

V'alue  (au) 

Notes 

xminl 

0.2 

minimum  X  grid  boundary 

0.2 

minimum  Y  grid  boundary 

2000.0 

mass  of  atom  a 

mb 

mass  of  atom  b 

me 

mass  of  atom  c 

X  grid  resolution 

deltax2 

0.05 

Y  grid  resolution 

nstep 

800 

number  of  time  steps 

sftent 

800 

number  of  time  steps  to  stop  the  propagation  at 

deltat 

1.0 

time  step 

ioffl 

1.4 

initial  X  offset  (0  for  product  wave  packet) 

iofF2 

0.0 

initial  Y  offset  (1.4  for  product  wave  packet) 

imoml 

-1.5 

initial  Px  offset  (0  for  product  wave  packet) 

imom2 

0.0 

initial  Py  offset  (-1.5  for  product  wave  packet) 

sprdl 

0.3 

initial  relative  wave  packet  spread  in  the  X  direction 

sprd2 

0.3 

initial  relative  wave  packet  spread  in  the  Y  direction 

taul 

-800.0 

temporal  offset  for  channel  1  initial  wave  packet 

tau2 

-800.0 

temporal  offset  for  channel  2  initial  wave  packet 

evch 

1 

asymptotic  eigenstate  number 

0  2  4  6  8  10  12 

X  (au) 


Figure  5.3  Compact  Mpller  states  and  absorbing  boundary  condition  placement 
for  the  two  coupled  Morse  oscillator  LLL  mass  configuration.  The  two 
parallel  contours  represent  the  0.1  au  contours  of  the  potential  energy 
surface.  The  dot  represents  the  bottom  of  the  well  in  the  interaction 
region.  The  two  closed  contours  are  the  0.1  au  contours  of  the  absolute 
value  of  the  Mpller  states.  The  shaded  region  represents  the  area  in 
which  the  absorbing  boundary  conditions  are  non-zero. 
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Absolute  Value  of  the  Wave  Packet 


Figure  5.4  Surface  plot  of  the  absolute  value  of  the  evolving  LLL  reactant  Mpller 
state  at  f  =  +800  au.  The  reactant  channel  is  parallel  to  the  X  direction 
and  the  product  channel  is  parallel  to  the  Y  direction.  Almost  all  the 
evolving  Moller  state  has  exited  the  interaction  region  in  the  product 
channel. 
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Absolute  Value  of  the  Correlation  Function 


0.14  1 


Relative  Translational  Energy  (eV) 

Figure  5.6  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  T  T  .1 ,  mass 
configuration  for  the  reaction  A  +  BC{v  =  0)  AB{y=  0) +C . 
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Figure  5.7  Absoute  value  of  the  correlation  function  for  the  two  coupled  Morse  oscil¬ 
lator  LLL  mass  configuration  for  the  reaction  A+BC{v  =  0)—> 

AB(v’=0)-i-C.  In  this  calculation,  the  ABC  have  been  moved  closer  to 
the  interaction  region  and  the  grid  size  reduced. 


Figure  5.8  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  LLL  mass 
configuration  for  the  reaction  A +SC(v  =  0) AB(v’=  0)+C.  In  this 

calculation,  the  ABC  have  been  moved  closer  to  the  interaction  region  and 
the  grid  size  reduced. 
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Figure  5.9  Absolute  value  of  the  correlation  function  for  the  two  coupled  Morse 
oscillator  LLL  mass  configuration  for  the  reaction  A  +  BC  (i/  =  0) 
AB  {v'  =  0)  +  C  where  absorbing  boundary  conditions  were  not  used. 

clearly  evident  in  the  probability  of  reaction  computed  from  the  correlation  function 
in  Figure  5.9  and  shown  in  Figure  5.10  and  is  likely  to  be  real.  In  Figure  5.10,  the 
probability  of  reaction  at  low  energies  and  the  overlying  oscillations  are  due  to  pre¬ 
maturely  terminating  the  computation  of  the  correlation  function  in  order  to  avoid 
edge  of  grid  reflection. 

5.2.3  MLM  S-matrix  dements.  Figure  5.11  illustrates  that  the  MLM 

Mpller  states  are  well  localized  in  the  interaction  region.  A  surface  plot  of  the 
evolving  reactant  Mpller  state,  shown  in  Figure  5.12,  illustrates  how  the  evolving 
Mpller  state  bifurcates  into  both  channels.  The  correlation  function,  shown  in  Fig¬ 
ure  5.13,  rapidly  approaches  zero  in  both  the  negative  and  positive  time  directions. 
However,  the  probability  of  reaction  shown  in  Figure  5.14  again  exhibits  an  overly¬ 
ing  oscillation  as  in  the  LLL  mass  conflguration  shown  in  Figure  5.6.  In  order  to 
determine  if  the  resonances  were  real  or  artifacts  of  absorbing  boundary  condition 
reflection  error,  the  correlation  function  was  computed  using  a  new  set  of  absorbing 
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Figure  5.10  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  LLL  mass 
configuration  for  the  reaction  A  +  BC  (i/  =  0)  — >■  AB  {u'  =  0)  +  C 
where  the  correlation  function  is  computed  without  the  use  of  absorb¬ 
ing  boundary  conditions. 

boundary  conditions,  shown  in  Figure  5.15,  and  without  using  absorbing  boundary 
conditions,  shown  in  Figure  5.16.  The  new  set  of  absorbing  boundary  conditions  are 
constructed  using  Jacobi  coordinates  appropriate  to  each  channel  and  then  trans¬ 
forming  them  to  bond  coordinates.  For  example,  in  the  A  -1-  BC  channel,  the  area  of 
the  channel  in  which  to  apply  absorbing  boundary  conditions  is  selected.  Then,  at 
each  point  in  this  area,  the  (Jf,  Y)  coordinates  are  transformed  according  to  Equation 
2.4  into  the  Jacobi  coordinates  {ra,Ra)-  Next,  the  absorbing  boundary  conditions 
are  computed  using  Equation  2.8  except  that  now  the  potential  is  given  in  terms  of 
Va  (ro,  Ra)  instead  of  K  (^)-  The  placement  of  the  absorbing  boundary  conditions 
is  still  subject  to  the  restriction  that  they  are  non-zero  only  where  the  initial  Mpller 
states  are  zero.  Figure  5.17  illustrates  the  probability  of  reaction  for  the  three  sets 
of  absorbing  boundary  conditions.  In  each  case,  the  oscillations  change  but  the  posi¬ 
tion  of  the  resonances  do  not.  The  resonances  and  the  underlying  shape  of  the  curve 
are  probably  real  wdth  the  overlying  oscillations  due  to  either  absorbing  boundary 
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Figure  5.11  Compact  Mpller  states  and  absorbing  boundary  condition  placement 
for  the  two  coupled  Morse  oscillator  MLM  mass  configuration.  The  two 
parallel  contours  represent  the  0.1  au  contours  of  the  potential  energy 
surface.  The  dot  represents  the  bottom  of  the  well  in  the  interaction 
region.  The  two  closed  contours  are  the  0.1  au  contours  of  the  absolute 
value  of  the  Mpller  states.  The  shaded  region  represents  the  area  in 
which  the  absorbing  boundary  conditions  are  non-zero. 
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Absolute  Value  of  the  Wave  Function 


Figure  5.12  Surface  plot  of  the  absolute  value  of  the  evolving  MLM  reactant  Mpller 
state  at  t  =  +1000  au.  The  reactant  channel  is  parallel  to  the  X 
direction  and  the  product  channel  is  parallel  to  the  Y  direction.  The 
evolving  Moller  state  has  bifurcated  into  both  channels  while  a  quasi¬ 
bound  state  is  evident  in  the  interaction  region. 


5-15 


Figure  5.13  Absolute  value  of  the  correlation  function  for  the  two  coupled  Morse 
oscillator  MLM  mass  configuration  for  the  reaction 

A+BC{v  =  0)  AB(v’=0)+C. 


Figure  5.14  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  MLM  mass 
configuration  for  the  reaction  A  +BC{y  =  0)  ^  AB{y=  0)+C. 
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Figure  5.15  Absoute  value  of  the  correlation  function  for  the  two  coupled  Morse  oscil¬ 
lator  MLM  mass  configuration  for  the  reaction  A -(-BC(v  =  0)  — > 
AB(v’=  0)  -t-C .  In  this  calculation,  the  ABC  are  first  constructed  in  Jacobi 
coordinates  and  then  transformed  to  bond  coordinates. 


Figure  5.16  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  MLM  mass 
configuration  for  the  reaction  A -f  5C(v  =  0)  — >  AB(v’=  0)-(-C.  In  this 

calculation,  the  ABC  are  first  constructed  in  Jacobi  coordinates  and  then 
transformed  to  bond  coordinates. 
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condition  reflection  error,  for  the  solid  and  dashed  lines,  or  are  due  to  prematurely 
terminating  the  computation  of  the  correlation  function  in  order  to  avoid  edge  of 
grid  reflection  for  the  dash-dot  line. 

5.2.4  HLH  S-matrix  elements.  Figure  5.18  illustrates  that  the  HLH 

Moller  states  are  well  localized  in  the  interaction  region.  In  the  correlation  function, 
shown  in  Figure  5.19,  the  onset  of  absorbing  boundary  condition  reflection  error  is 
quite  clearly  visible  at  approximately  t  =  -i-6000  au.  In  order  to  test  that  the  oscil¬ 
lations  in  Figure  5.19  are  due  to  absorbing  boundary  condition  reflection,  another 
calculation  was  performed  using  the  same  absorbing  boundary  conditions  placed  fur¬ 
ther  from  the  interaction  region  on  a  larger  grid.  Figure  5.20  illustrates  that  as  the 
absorbing  boundary  conditions  are  moved  further  from  the  interaction  region,  the 
onset  of  oscillations  in  the  correlation  function  do  occur  later  in  time  are  are  due  to 
absorbing  boundary  condition  reflection  error.  The  probability  of  reaction,  shown 
in  Figure  5.21  using  the  correlation  function  shown  in  Figure  5.19,  e.xhibits  a  reso¬ 
nance  at  0.34  eV  and  an  overlying  oscillation.  As  in  both  the  LLL  and  MLM  mass 
configurations,  another  computation  of  the  correlation  function  without  the  use  of 
absorbing  boundary  conditions  was  performed  in  order  to  determine  which  features 
in  the  probability  of  reaction  in  Figure  5.21  are  most  likely  real  and  which  are  due 
to  reflection  error.  Figure  .5.22  illustrates  the  resulting  probability  of  reaction  for  the 
three  calculations.  Figure  5.23  illustrates  the  correlation  function  where  the  compu¬ 
tation  does  not  use  absorbing  boundary  conditions.  The  oscillations  at  t  =  -1-6000 
au  are  absent.  However,  in  order  to  avoid  edge  of  grid  reflections,  the  computation 
is  terminated  prematurely  which  will  lead  to  spurious  oscillations  in  the  resulting  S- 
matrix  elements.  In  each  calculation,  the  resonance  in  the  probability  of  reaction  at 
0.34  eV  remains  and  is  most  likely  real  as  well  as  the  underlying  shape  of  the  curve. 
However,  the  remaining  oscillations  change  in  both  amplitude  and  frequency  are  are 
due  to  errors  introduced  by  reflection  from  the  absorbing  boundary  conditions  (solid 
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Figure  5.18  Compact  Mpller  states  and  absorbing  boundary  condition  placement 
for  the  two  coupled  Morse  oscillator  HLH  mass  configuration.  The  two 
parallel  contours  represent  the  0.1  au  contours  of  the  potential  energy 
surface.  The  dot  represents  the  bottom  of  the  well  in  the  interaction 
region.  The  two  closed  contours  are  the  0.1  au  contours  of  the  absolute 
value  of  the  Mpller  states.  The  shaded  region  represents  the  area  in 
which  the  absorbing  boundary  conditions  are  non-zero. 
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Figure  5.20  Absolute  value  of  the  correlation  function  for  the  two  coupled  Morse 
oscillator  HLH  mass  configuration  for  the  reaction 
A  +  BC{v  =  0)  ->  AB{v'=  0)+C  on  two  different  grid  sizes. 


Relative  Translational  Energy  (eV) 


Figure  5.21  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  HLH  mass 
configuration  for  the  reaction  A  +  BC(v  =  0) -4  AB(v’=  0)+C.  The 
correlation  function  used  to  compute  this  probability  is  shown  in  Figure 
5.19. 
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Figure  5.23  Absolute  value  of  the  correlation  function  for  the  two  coupled  Morse 
oscillator  HLH  mass  configuration  for  the  reaction  A  +  BC  (i/  =  0)  -» 
AB  {u'  =  0)  +  C  where  absorbing  boundary  conditions  were  not  used. 

and  dashed  lines)  or  premature  termination  of  the  correlation  function  calculation 
(dash-dot  line). 

5.3  S-matrix  elements  for  absorbing  boundary  condition  reflection 

In  order  to  more  fully  investigate  the  problem  of  absorbing  boundary 
condition  reflection  error,  S-matrix  elements  were  computed  for  absorbing  boundary 
condition  reflection  using  a  simple  trough  potential  energy  surface  for  the  LLL,  MLM, 
and  HLH  mass  configurations. 

5.3.1  Trough  potential  energy  surface.  The  potential  energy  surface  is 

simply  a  trough  parallel  to  the  X  direction  and  a  Morse  oscillator  in  the  Y  direction 
given  by 

Vtrough  (^,  y)  =  De  {I-  CXp  [-«  {Y  -  Te)]}^  ,  (5.4) 

where  Vtrough  has  the  same  form  as  the  asymptotic  limit  of  the  channels  in  the  full  two 
dimensional  calculation.  In  this  calculation,  the  scattering  potential  (which  defines 
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the  interaction  region)  is  chosen  to  be  the  absorbing  boundary  conditions  used  in 
the  full  two  dimensional  calculation.  Because,  in  the  channel  packet  method,  we 
are  free  to  choose  the  coordinate  offset  and  zero  of  time,  the  initial  wave  packets 
and  final  Moller  states  can  be  the  same  by  placing  the  coordinate  offset  far  from 
the  absorbing  boundary  conditions.  Figure  5.24  illustrates  these  points.  The  shaded 
region  represents  the  area  in  which  the  absorbing  boundary  conditions  are  non-zero. 
In  this  calculation,  the  absorbing  boundary  conditions  are  the  scattering  potential 
and  thus  also  represent  the  interaction  region.  The  two  parallel  lines  represent  the 
0.1  au  contours  of  the  trough  potential  given  in  Equation  5.4.  The  closed  contours  at 
the  right  end  of  the  potential  represent  both  the  initial  wave  packet  and  final  Moller 
state.  They  are  the  same  since  the  initial  wave  packet  is  placed  in  the  asymptotic 
limit  and  under  propagation  by  the  Mpller  Operators  given  in  Equation  2.10,  the 
asymptotic  and  full  scattering  Hamiltonians  are  the  same  thus  the  effect  of  the  Moller 
Operator  is  multiplication  by  1. 

5.3.2  S-matrix  elements.  Figures  5.25,  5.26  and  5.27  illustrate  the 

probability  of  reflection  for  the  reaction  .4  -|-  BC  {v  =  j)  ^  A  BC  {u'  =  j) ,  j  = 
0, 1,2  for  the  LLL,  MLM  and  HLH  mass  configurations.  Comparison  of  the  overall 
magnitudes  of  the  probability  of  reflection  for  the  three  mass  configurations  and, 
within  a  given  mass  configuration,  the  three  internal  vibrational  eigenstates,  indi¬ 
cates  that  wave  packets  in  the  MLM  mass  configuration  are  less  likely  to  undergo 
absorbing  boundary  condition  reflection  than  the  LLL  and  HLH  mass  configura¬ 
tions.  For  some  energies,  the  HLH  or  LLL  mass  configurations  may  have  a  lower 
probability  of  reflection,  overall,  they  are  more  likely  to  reflect  than  the  MLM  mass 
configuration.  In  the  MLM  and  HLH  mass  configurations  wave  packets  constructed 
from  the  u  =  0  internal  vibrational  eigenstate  are  more  likely  to  undergo  absorbing 
boundary  condition  reflection  than  wave  packets  constructed  from  the  u  =  1  and 
u  =  2  internal  vibrational  eigenstates.  The  probability  of  reflection  for  t/  ^  u  is 
essentially  zero. 
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Summary 

The  initial  purpose  of  this  part  of  the  study  was  to  investigate,  using  the 
channel  packet  method  with  absorbing  boundary  conditions,  the  effects  of  potential 
energy  surface  w’ells  on  the  probability  of  reaction.  However,  absorbing  boundary 
condition  reflections  introduce  non-trivial  errors  into  the  correlation  function  with 
resulting  errors  in  the  S-matrix  elements.  Currently,  the  unambiguous  signatures 
of  absorbing  boundary  condition  reflection  error  are  S-matrix  elements  that  fail  to 
converge,  probabilities  of  reaction  significantly  greater  than  one  or  oscillations  in  the 
probability  of  reaction  which  change  with  changes  in  the  absorbing  boundary  con¬ 
ditions.  Some  insight  into  which  features  in  the  probability  of  reaction  are  real  and 
which  are  due  to  absorbing  boundary  condition  reflections  can  by  obtained  by  com¬ 
puting  the  correlation  function  without  absorbing  boundary  conditions.  However, 
these  comparisons  only  yield  insight  into  where  there  might  be  sharp  resonances  and 
the  overall  trends  in  the  probability  of  reaction. 
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VI.  Two  dimensional  reactive  quantum  scattering  for  a  model 
system  of  two  coupled  Morse  oscillators  in  a  light-heavy-light 

configuration 

Quantum  reactive  scattering  results  are  presented  for  a  model  system  of 
two  coupled  Morse  oscillators  in  a  light-heavy-light  (LHL)  mass  configuration.  The 
conclusions  in  this  chapter  are  based  on  observed  features  of  the  S-matrix  elements 
and  are  primarily  qualitative  in  nature. 

6.1  S-matrix  dements 

The  possibility  of  significant  absorbing  boundary  condition  reflection  error 
for  the  LHL  mass  configuration  needs  to  be  eliminated  before  the  effects  of  well  depth 
and  kinetic  coupling  may  be  explored.  For  the  LHL  mass  configuration,  the  light 
masses  are  2000  au  and  the  heavy  mass  is  10,000  au.  Using  the  same  model  potential 
energy  surface  as  in  Equation  5.1  with  Dg  =  0.1  au  and  a  —  5.0,  the  asymptotic 
limit  of  the  Morse  oscillator  has  4  bound  states  where  the  i/  =  0,1,2  vibrational 
eigenstate  energies  are  the  same  as  for  the  HLH  mass  configuration  in  Table  5.2. 
The  energies  are  the  same  since  for  both  the  HLH  and  LHL  mass  configurations  the 
diatom  mass  configuration  consists  of  a  light  atom  and  a  heavy  atom.  Parameters  for 
the  Mpller  state  calculation  are  given  in  Table  6.1  and  the  procedure  for  calculating 
the  Mpller  states  is  the  same  as  in  Section  5.2.1.  The  probability  of  reaction  shown 
in  Figure  6.1  exhibits  three  sharp  resonances.  However,  in  light  of  the  results  in 
Chapter  V,  where  absorbing  boundary  condition  reflection  introduced  non-trivial 
errors  into  the  calculation,  two  further  computations  of  the  correlation  function  were 
performed  in  order  to  investigate  the  effects,  if  any,  of  absorbing  boundary  condition 
reflection  error.  The  first  computation  involved  a  larger  grid  where  the  absorbing 
boundary  conditions  were  moved  further  from  the  interaction  region.  Figure  6.2 
illustrates  the  results  compared  to  the  correlation  function  computation  on  a  smaller 
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Table  6.1  Initial  M0ller  state  calculation  parameters  for  the  LHL  mass 
configuration. 


Parameter 

Value  (au) 

Notes 

xminl 

0.2 

minimum  X  grid  boundary 

xmin2 

0.2 

minimum  Y  grid  boundary 

ma 

2000.0 

mass  of  atom  a 

mb 

10000.0 

mass  of  atom  b 

me 

2000.0 

mass  of  atom  c 

deltaxl 

0.05 

X  grid  resolution 

deltax2 

Y  grid  resolution 

nstep 

number  of  time  steps 

sftent 

1200 

number  of  time  steps  to  stop  the  propagation  at 

deltat 

1.0 

time  step 

ioffl 

1.4 

initial  X  offset  (0  for  product  wave  packet) 

ioff2 

0.0 

initial  Y  offset  (1.4  for  product  wave  packet) 

imoml 

-1.4 

initial  Px  offset  (0  for  product  wave  packet) 

imom2 

0.0 

initial  Py  offset  (-1.4  for  product  wave  packet) 

sprdl 

0.3 

initial  relative  wave  packet  spread  in  the  X  direction 

sprd2 

0.3 

initial  relative  wave  packet  spread  in  the  Y  direction 

temporal  offset  for  channel  1  initial  wave  packet 

temporal  offset  for  channel  2  initial  wave  packet 

evch 

1 

asymptotic  eigenstate  number 
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Figure  6.2  Absolute  value  of  the  correlation  function  for  the  two  coupled  Morse 
oscillator  LHL  mass  configuration  for  the  reaction 
A  +  BC(v  =  0)  — >  AB(v’=  0)+C  for  two  different  grid  sizes.  For  the 
larger  grid  size,  the  absorbing  boundary  conditions  are  placed  further  away 
from  the  interaction  region. 
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Figure  6.3  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  LHL  mass 
configuration  for  the  reaction  A-\-  BC  (0)  — >  AB  (0)  +  C  for  two  different 
grid  sizes.  For  the  larger  grid  size,  the  absorbing  boundary  conditions 
are  placed  further  away  from  the  interaction  region. 


grid.  The  differences  between  the  two  functions  are  small  and  the  resulting  S- 
matrix  elements,  shown  in  Figure  6.3,  are  essentially  the  same  as  in  Figure  6.1.  The 
second  computation,  shown  in  Figure  6.4,  involved  using  the  large  grid  without  the 
application  of  absorbing  boundary  conditions.  The  computation  of  the  correlation 
function  without  absorbing  boundary  conditions  introduces  oscillations  into  the  S- 
matrix  elements  shown  in  Figure  6.5.  These  oscillations  are  due  to  terminating  the 
correlation  function  calculation  before  the  correlation  function  has  reached  zero  in 
order  to  avoid  the  onset  of  edge  of  grid  reflection.  The  comparison  of  all  three 
computations  demonstrates  that  absorbing  boundary  condition  reflection  introduces 
negligible  error  into  the  resulting  S— matrix  elements  for  the  LHL  mass  configuration. 

The  original  intent  of  combining  absorbing  boundary  conditions  with  the  chan¬ 
nel  packet  method  was  to  improve  the  efficiency  of  the  computation  of  S-matrix 
elements.  Figure  6.5,  where  the  correlation  function  is  computed  without  the  use  of 
absorbing  boundary  conditions,  illustrates  that  on  a  512x512  grid  without  absorbing 
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Figure  6.4  Absolute  value  of  the  correlation  function  for  the  two  coupled  Morse 
oscillator  LHL  mass  configuration  for  the  reaction 
A  +  BC{v  =  0)  — >  AB[v’=  0)+C  without  the  use  of  absorbing  boundary 

conditions. 


Relative  Translational  Energy  (eV) 


Figure  6.5  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  LHL  mass 
configuration  for  the  reaction  A+BC{v  =  0) —>  AB[v’=0)+C  without 
the  use  of  absorbing  boundary  conditions. 
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Figure  6.6  Absolute  value  of  the  correlation  function  for  the  two  coupled  Morse 
oscillator  LHL  mass  configuration  for  the  reaction  A  +  BC  (0)  — >• 
AS (0)  +  C  where  the  grid  size  has  been  reduced  to  128x128  and  the 
absorbing  boundary  conditions  moved  closer  to  the  interaction  region. 

boundary  conditions,  good  S-matrix  elements  can  be  obtained.  The  calculation  of 
the  correlation  function  in  Figure  6.4,  on  a  MIPS  RIOOOO  Silicon  Graphics  work¬ 
station  took  approximately  55,000  seconds  of  CPU  time.  However,  using  absorbing 
boundary  conditions  and  a  small  128x128  grid,  the  calculation  of  the  correlation 
function,  shown  in  Figure  6.6,  took  only  1900  seconds  of  CPU  time  for  a  factor  of 
29  savings  in  computation  time.  As  shown  in  Figure  6.7,  the  resulting  S-matrix  el¬ 
ements  computed  from  the  correlation  function  in  Figure  6.6  are  in  good  agreement 
w'ith  Figures  6.1  and  6.3.  Again,  the  combination  of  absorbing  boundary  conditions 
with  the  channel  packet  method  results  in  significant  time  savings  as  opposed  to  the 
channel  packet  method  alone. 

6.2  Resonances  and  well  depth 

Since  the  LHL  mass  configuration  is  effectively  free  of  absorbing  boundary 
condition  reflection  error,  the  influence  of  well  depth  on  S-matrix  elements  can  be 
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Figure  6.7  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  LHL  mass 
configuration  for  the  reaction  A  +  BC  (0)  — >•  AB  (0)  +  C  where  the 
grid  size  has  been  reduced  to  128x128  and  the  absorbing  boundary 
conditions  moved  closer  to  the  interaction  region. 

investigated.  The  well  depth  in  the  model  potential  energy  surface  is  easily  adjusted 
by  changing  the  dissociation  energy  of  the  Morse  oscillators,  denoted  by  in  Equa¬ 
tion  5.1.  In  this  investigation,  S-matrix  elements  are  computed  for  values  of  £>«  from 
0.1  au  to  0.26  au  in  steps  of  0.01  au.  The  step  size  is  purely  arbitrary.  For  each  value 
of  De,  the  anharmonicity  a  was  changed  so  that  the  spring  constant  of  the  harmonic 
limit  of  the  Morse  oscillator,  given  by  K  =  2Deoi^,  remained  the  same.  For  each 
potential  energy  surface,  the  initial  grid  and  wave  packet  parameters  are  the  same 
as  in  Table  6.1.  Figure  6.8  illustrates  the  probability  of  reaction  for  four  values  of 
De-  Figures  illustrating  the  probability  of  reaction  for  values  of  Dg  from  0.1  au  to 
0.26  au  in  steps  of  0.01  au  are  in  Appendix  D. 

At  first  glance,  a  higher  dissociation  energy,  i.e.  a  deeper  well,  leads  to  more 
resonances  with  doublets  and  then  triplets  appearing.  By  comparing  S-matrix  el¬ 
ements  computed  using  different  values  of  Dg.  resonances  that  initially  appear  at 
high  relative  translational  energy  move  to  lower  relative  translation  energy  as  the 
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Probability  of  Reaction 


Figure  6.8  Probability  of  reaction  as  a  function  of  dissociation  energy  for  the  two 
coupled  Morse  oscillator  LHL  mass  configuration  for  the  reaction  A  + 
BC  (0)  —>  AB  (0)  +  C.  The  dissociation  energy  is  equal  to  0.10  au,  0.15 
au,  0.20  au  and  0.26  au. 


dissociation  energy  increases.  The  peak  positions  of  several  of  the  resonances  with 
respect  to  potential  energy  surface  dissociation  energy  is  shown  in  Figure  6.9.  The 
initial  position  appears,  qualitatively,  to  be  governed  by  an  envelope.  The  form  of 
the  envelope  can  be  fit  by  a  function  of  the  form 


(6.1) 


where  y  is  the  resonance  position  and  x  is  the  dissociation  energy.  This  envelope  as 
well  as  the  points  used  to  derive  it  are  shown  in  Figure  6.10. 

In  Figure  6.11,  the  position  of  a  resonance  with  respect  to  dissociation  energy 
is  connected  with  a  curve  for  a  number  of  resonances.  The  connections  between  the 
points  in  Figure  6.9  which  are  not  included  in  Figure  6.11  are  more  difficult  to  make 
due  to  the  bifurcation  and  trifurcation  of  the  underlying  resonances.  Each  of  the 
curves  in  Figure  6.11  has  the  same  general  shape,  wdth  the  position  of  a  resonance 
changing  slowly  at  first  and  then  more  rapidly  with  increasing  dissociation  energy. 
Some  resonances  also  appear  to  bifurcate  with  increasing  well  depth.  In  Figure  6.25, 
the  resonance  represented  by  the  dashed  line  was  chosen  for  further  investigation 
due  to  the  bifurcation  at  a  dissociation  energy  of  0.19  au.  S-matri.x  elements  were 
computed  for  =  0.18  to  0.20  au  in  steps  of  0.002  au.  Figure  6.12  illustrates 
the  position  of  the  bifurcating  resonance  as  a  function  of  dissociation  energy.  The 
bifurcation  occurs  between  0.184  and  0.186  au  and  the  subsequent  behavior  of  the 
resonance  position  rapidly  approaches  the  general  behavior  of  the  other  resonances 
as  seen  in  Figure  6.9. 

For  the  one  dimensional  square  well,  the  resonances  in  the  probability  of  trans¬ 
mission  can  be  related  to  virtual  states  of  the  infinitely  deep  square  well.^®’'*'  In  order 
to  investigate  the  possibility  that  such  a  relationship  might  exist  for  the  two  coupled 
Morse  oscillators,  the  eigenstates  of  the  harmonic  limit  need  to  be  computed.  This 
can  be  done  analytically  by  constructing  the  force,  F,  and  mass,  M,  matrices  then 
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Figure  6.9  Position  of  S-matrix  element  resolances  as  a  function  of  Morse  oscillator  dissociation  energy. 


Figure  6.10  Qualitative  envelope  governing  the  initial  position  of  a  resonance  in 
the  probability  of  reaction. 


diagonalizing  the  matrix  (M"*F).“  The  diagonal  elements  of  the  resulting  matrix 
are  the  squares  of  the  frequencies  of  the  normal  modes  of  the  system.  For  the  two 
coupled  Morse  oscillator  potential,  the  force  matrix  is  given  by 


F  = 


K 

-K 

0 


-K 

2K 

-K 


0  \ 

-K 

K  J 


(6.2) 


where  K  is  the  force  constant  of  the  harmonic  limit  of  the  Morse  oscillator  and  the 
mass  matrix  is  diagonal  with  rria  =  rric  and  given  by 
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Diagonalizing  the  matrix  (M  ^F)  yields 


/^o 

0 

0  \ 

/o 

0 

0  \ 

0 

0 

0 

iL 

771a 

0 

0 

0 

Uj]} 

u 

0 
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rria^b 

The  normal  mode  associated  with  wo  =  0  is  translation.  The  normal  mode  associated 
with  uJi  is  the  symmetric  stretch  and  the  normal  mode  associated  with  ^>2  is  the 
asymmetric  stretch.  The  energy  eigenvalues  of  the  harmonic  limit  of  the  two  coupled 
Morse  oscillator  potential  are  given  by 


E  (ni,n2)  =  ^^1 


flLC2 


(6.5) 


where  rii  and  n2  are  the  quantum  numbers  associated  with  the  eigenstates  (normal 
modes)  whose  frequencies  are  oji  and  uj2-  «i  and  n2  are  integers  and  range  from 
zero  to  infinity.  For  the  LHL  mass  configuration,  the  energy  eigenvalues  associated 
with  the  frequencies  are  =  1.36  eV  and  hu!2  =  1.61  eV.  The  20  lowest  energy 
eigenvalues  are  given  in  Table  6.2  with  their  associated  quantum  numbers.  The 
spectrum  of  the  energy  eigenvalues  is  given  in  Figure  6.13.  From  the  spectrum 
below  8  eV,  the  eigenstates  of  the  harmonic  limit  can  be  grouped  into  singlets, 
doublets  and  triplets.  However,  the  resonances  in  the  S-matrix  elements  shown  in 
Figures  6.1  and  6.8  do  not  line  up  with  the  spectrum  shown  in  Figure  6.13,  and 
there  is  no  clear  analytic  connection  between  the  resonances  and  the  eigenstates  of 
the  harmonic  limit.  Graphically,  at  low  values  of  dissociation  energy,  the  resonances 
can  be  grouped  as  singlets  and  doublets.  At  higher  dissociation  energy,  resonances 
can  also  be  grouped  into  triplets.  It  is  tempting  to  attribute  the  resonant  multiplets 
to  virtual  states  of  the  harmonic  limit  of  the  two  coupled  Morse  oscillators.  However, 
only  qualitative  similarities  are  observed  and  making  such  a  connection  between  the 
virtual  states  of  the  harmonic  limit  and  the  resonances  in  Figures  6.1  and  6.8  is 
premature. 
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6.3  Resonances  and  kinetic  coupling 

In  the  bond  coordinate  representation  of  the  full  scattering  Hamiltonian, 
the  kinetic  energy  operators  in  Equation  5.2  are  coupled  via 


PxPy 

mi, 


(6.6) 


w'here  Pj  is  the  conjugate  momentum  operator  in  the  j  direction  and  mi,  is  the  mass 
of  atom  B.  In  this  analysis,  the  minus  sign  in  front  of  the  coupling  term  is  replaced 
bv  a  parameter  which  is  varied  from  -1  to  0  in  steps  of  0.1  where  the  step  size  is 
arbitrary.  The  potential  energy  surface  with  a  dissociation  energy  of  0.20  au  was 
chosen  for  the  analysis  for  the  number  of,  and  complex  nature  of,  the  resonances. 
Figures  6.14  and  6.15  illustrate  the  effect  of  kinetic  coupling  on  S-matrix  resonances. 
The  complete  set  of  S-matrix  elements  for  kinetic  energy  coupling  constants  of-1  to 
0  in  steps  of  0.1  are  contained  in  Appendix  D.  As  the  kinetic  coupling  is  reduced,  the 
resonances  labeled  1  in  Figure  6.14,  simply  decrease  in  magnitude  with  decreasing 
coupling.  The  resonances  labeled  2  both  decrease  in  magnitude  and  coalesce  into 
a  single  resonance.  .\s  expected,  the  probability  of  reaction  where  the  two  Morse 
oscillators  are  completely  kinetically  decoupled  is  zero  as  shown  in  Figure  6.15.“ 


6-4  Absorbing  boundary  condition  reflection  S-matrix  elements 

Similar  to  the  analysis  in  Section  5.3,  S-matrix  elements  were  computed 
for  absorbing  boundary  condition  reflection  for  the  reaction  A  +  BC  {i/  =  j)  A  A 
BC  {u'  =  j) ,  j  =  0, 1, 2  for  the  LHL  mass  configuration  and  are  shown  in  Figure  6.16. 
The  probability  of  reflection  for  i/'  u  \s  essentially  zero.  As  seen  in  the  MLM 
and  HLH  mass  configurations,  wave  packets  constructed  from  the  u  =  0  internal 
vibrational  eigenstate  are  more  likely  to  reflect  than  those  constructed  from  the 
1/  =  1  and  u  =  2  eigenstates.  However,  the  magnitude  of  the  probability  of  reflection 

“The  labels  1  and  2  are  purely  arbitrary. 
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is  not  drastically  reduced  from  the  probability  of  reflection  for  the  MLM  or  HLH 
cases. 

The  reason  why  the  LHL  mass  conflguration  does  not  suffer  from  significant 
absorbing  boundary  condition  reflection  error  is  not  attributable  solely  to  the  magni¬ 
tude  of  absorbing  boundary  condition  reflection.  In  order  to  investigate  the  reasons 
why  absorbing  boundary  condition  reflection  error  is  significant  in  the  LLL,  MLM 
and  HLH  mass  configurations  and  not  the  LHL  mass  configurations,  consider  the 
correlation  function  given  by 

Cy,  (()  =  (^z'lexp  m>-  («•') 

The  derivation  of  the  formula  for  S-matrix  elements  in  Appendix  A  and  Chapter 
II  uses  a  Hamiltonian  that  does  not  contain  absorbing  boundary  conditions.  With 
the  application  of  absorbing  boundary  conditions,  the  formula  for  the  correlation 
function  changes, 

Cy,(()  =  {«l’|exp{^  (h  -  (6-8) 

to  include  K  the  absorbing  boundary  condition  potential.  The  wave  packets  |^)  can 
be  expanded  in  terms  of  expansion  coefficients  and  the  eigenbasis  of  H , 

/'  +  00 

l«i>=/  (6.9) 

J—oo 

Expanding  the  bra-ket  in  Equation  6.8  and  using  Equation  6.9,  the  result  is 

C-y'-y  {t)  =  f  dky  [  dk^rjl  (fey)  ry+  {k^) 

J —oo  J—oo 

x(/i:y,7'-|exp|-^  (^H  -  iVa)^\k-yr/+)-  (6-10) 
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VVe  believe  that  the  mixing  of  eigenstates  of  H  by  K  in  Equation  6.10  governs 
the  absorbing  boundary  condition  reflection  error.  However,  the  complex  nature  of 
the  integrals  in  Equation  6.10  does  not  lend  itself  to  deriving  an  absorbing  potential 
which  minimizes  the  error  introduced  by  mixing  the  eigenstates  of  H  and  a  solution 
to  the  problem  of  absorbing  boundary  condition  reflection  error  is  unknown  at  this 
point. 

6.5  Summary 

The  light-heavy-light  mass  configuration  of  the  two  coupled  Morse  oscilla¬ 
tors  is  qualitatively  unaffected  by  absorbing  boundary  condition  error.  This  can  be 
seen  in  that  changing  the  position  of  the  absorbing  boundary  conditions  does  not  af¬ 
fect  the  probability  of  reaction.  This  is  in  contrast  to  the  results  in  Chapter  V  where 
the  position  of  the  absorbing  boundary  conditions  greatly  influence  the  resulting 
probability  of  reaction.  One  of  the  primary  purposes  of  this  study  was  to  improve 
the  efficiency  of  the  computation  of  S-matrix  elements.  For  the  light-heavy-light 
mass  configuration,  the  correlation  function  computation  without  absorbing  bound¬ 
ary  conditions,  which  yields  good  S-matrix  elements,  took  approximately  .55,000 
seconds  of  CPU  time.  Applying  absorbing  boundary  conditions,  reduced  the  grid 
size  by  a  factor  of  16  yielding  a  factor  of  29  reduction  in  the  computation  time  for 
the  correlation  function.  Again,  the  combination  of  absorbing  boundary  conditions 
with  the  channel  packet  method  results  in  significant  time  savings  as  opposed  to  the 
channel  packet  method  alone. 

The  effect  of  dissociation  energy  on  the  probability  of  reaction  for  the  two  cou¬ 
pled  Morse  oscillators  can  be  investigated  using  the  light-heavy-light  mass  configu¬ 
ration  since  it  is  not  affected  by  absorbing  boundary  condition  error.  Qualitatively, 
the  deeper  the  well  in  the  interaction  region,  i.e.  the  higher  the  dissociation  energy, 
the  probability  of  reaction  exhibits  more  resonances.  The  initial  placement,  energet¬ 
ically,  of  a  resonance  appears  to  be  governed  by  an  envelope  shown  in  Figure  6.10. 
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Overall,  the  position  of  a  resonance,  as  seen  in  Figure  6.11,  as  a  function  of  dissoci¬ 
ation  energy  exhibits  the  same  overall  shape.  Figure  6.12  illustrates  the  breaking  of 
a  degenerate  resonance  between  0.184  and  0.186  au. 

The  virtual  eigenstate  of  the  harmonic  limit  of  two  coupled  Morse  oscillators 
can  be  determined.  The  energy  spectrum  of  these  states  is  shown  in  Figure  6.1.3. 
However,  no  clear  analytic  connection  can  be  made  between  virtual  eigenstates  of  the 
harmonic  limit  of  the  two  coupled  Morse  oscillators  and  the  resonances  illustrated  in 
Figures  6.1  and  6.8.  The  question  of  influence  by  the  virtual  states  remains  open. 

The  effect  of  kinetic  coupling  is  also  investigated  using  the  light-heavy-light 
mass  configuration.  Figures  6.14  and  6.15  illustrate  the  effect  of  kinetic  coupling  on 
S-matrix  resonances.  .\s  the  kinetic  coupling  constant  is  reduced,  so  is  the  magnitude 
and  number  of  S-matrix  resonances.  When  the  two  Morse  oscillators  are  completely 
uncoupled,  the  coupling  constant  is  zero,  the  S-matrix  elements  are  zero.  This  is 
expected  since  without  any  coupling  there  is  no  energy  transfer  between  the  two 
oscillators.  For  a  tiny  coupling  constant,  -0.1,  there  is  only  a  slight  possibility  of 
reaction  at  a  relative  translation  energy  of  approximately  0.5  eV.  -A.n  analysis  of  S- 
matrix  resonances  divides  the  resonances  into  two  types.  One  type  simply  decreases 
in  magnitude  with  decreasing  coupling.  The  other  type  concerns  pairs  of  resonances 
that  coalesce  into  a  single  degenerate  resonance  before  their  magnitude  reaches  zero 
with  decreasing  coupling. 
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VII.  Conclusion 


The  calculation  of  S-matrix  elements  for  quantum  reactive  scattering  re¬ 
mains  computationally  intensive.  VVe  have  developed  a  new  method  that  dramati¬ 
cally  improves  efficiency  for  computing  S-matrix  elements  by  combining  the  channel 
packet  method  with  absorbing  boundary  conditions.  While  absorbing  boundary 
conditions  have  been  used  in  the  past  for  preventing  edge  of  grid  reflections,  never 
before  have  they  been  used  for  the  specific  purpose  of  reducing  the  size  of  the  grid 
used  to  perform  the  computations.  The  dramatic  improvement  in  computational 
efficiency  is  achieved  by  e.xploiting  the  manner  in  which  the  channel  packet  method 
computes  S-matrix  elements.  The  channel  packet  method  relies  on  computing  two 
Mpller  states:  one  representing  reactants  and  one  representing  products.  The  time 
dependent  correlation  function  between  the  two  Mpller  states  is  computed  and  used 
to  calculate  S-matrix  elements.  The  formulation  using  two  Mpller  states  is  exploited 
by  individually  propagating  the  Mpller  states  on  a  small  efficient  grid  by  using  ab¬ 
sorbing  boundary  conditions  to  attenuate  the  evolving  Mpller  states  as  they  exit  the 
interaction  region  of  the  potential.  The  attenuation  of  the  evolving  states  prevents 
their  reflecting  from  the  edges  of  the  much  smaller  grid.  The  Fourier  transform  of 
the  correlation  function  is  then  used  to  compute  S-matrix  elements. 

The  one  dimensional  square  well  has  an  analytic  solution  for  the  S-matrix 
elements.  S-matrix  elements  computed  using  the  combination  of  the  channel  packet 
method  with  absorbing  boundary  conditions  are  in  excellent  agreement  with  the 
analytic  solution.  Two  tests  of  convergence  are  used  to  demonstrate  agreement  with 
the  analytic  solution.  First,  it  is  important  that  a  method  not  only  converge  but  that 
it  converge  to  the  correct  answer.  In  Figure  3.6,  the  solution  computed  using  the 
channel  packet  method  with  absorbing  boundary  conditions  is  converging  towards 
the  analytic  solution.  Second,  the  order  of  convergence  provides  a  test  of  efficiency 
for  the  channel  packet  method  with  absorbing  boundary  conditions.  The  order  of 
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convergence  can  be  computed  by  determining  the  slope  of  a  line  in  a  log-log  plot 
where  the  changes  in  a  given  parameter,  for  example  the  spatial  grid  resolution, 
are  plotted  against  the  relative  error.  For  changes  in  the  spatial  grid,  the  channel 
packet  method  with  absorbing  boundary  conditions  converges  with  an  order  of  3. 
For  changes  in  the  temporal  grid,  the  order  is  4.  These  orders  of  convergence  are 
illustrated  in  Figures  3.8  and  3.9.  These  convergence  tests  validate  the  application  of 
absorbing  boundary  conditions  with  the  channel  packet  method  to  one  dimensional 
quantum  reactive  scattering. 

While  absorbing  boundary  conditions  are  extensively  used,  very  little  informa¬ 
tion  is  available  about  their  behavior  and  possible  influence  on  S-matrix  elements. 
Quantum  scattering  for  a  simple  one  dimensional  potential  consisting  of  a  Gaussian 
well  bounded  by  symmetric  Gaussian  barriers  is  used  to  investigate  the  behavior  of 
absorbing  boundary  conditions  when  applied  to  the  channel  packet  method.  The 
efficiency  of  the  combination  is  once  again  demonstrated  when  quasi-bound  states 
are  trapped  by  the  barriers  around  the  Gaussian  well.  Without  absorbing  boundary 
conditions,  the  grid  would  have  to  be  large  enough  to  support  the  evolving  Mpller 
state  until  the  quasi-bound  state  has  completely  exited  the  interaction  region.  For 
example,  when  the  barriers  are  0.0.5  au  high,  without  absorbing  boundary  conditions, 
the  spatial  grid  would  have  to  be  21  times  larger  in  order  to  accommodate  the  evolv¬ 
ing  wave  packet  until  the  quasi-bound  state  has  completely  exited  the  interaction 
region.  The  resulting  large  grid  is  computationally  prohibitive.  The  investigation 
into  the  effects  of  absorbing  boundary  conditions  begins  with  conditions  that  are  too 
steep.  Absorbing  boundary  conditions  that  are  too  steep  will  reflect  the  evolving 
Mpller  state  back  into  the  interaction  region  thereby  introducing  error  into  the  cor¬ 
relation  function.  The  resulting  S-matrix  elements  exhibit  near  periodic  oscillations 
and  probabilities  of  transmission  greater  than  one.  Similarly,  absorbing  boundary 
conditions  which  are  too  shallow  will  not  fully  attenuate  the  evolving  Mpller  state 
and  the  periodic  boundary  conditions  imposed  by  the  fast  Fourier  transforms  used 
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in  the  propagation  allow  the  Mpller  state  to  wrap  around  the  grid.  This  introduces 
errors  into  the  correlation  function  and  the  resulting  S-matrix  elements  again  exhibit 
nearly  periodic  oscillations  as  well  as  probabilities  of  transmission  greater  than  one. 
Finally,  the  sudden  onset  of  the  absorbing  boundary  conditions  can  introduce  small 
errors  into  the  correlation  function.  The  resulting  S-matrix  elements  exhibit  only 
minor  oscillations  about  the  correct  answer. 

Quantum  reactive  scattering  for  the  collinear  H  H2  reaction  has  been  ex¬ 
tensively  studied  using  a  variety  of  methods  and  provides  an  excellent  benchmark 
for  testing  the  combination  of  the  channel  packet  method  with  absorbing  boundary 
conditions.  The  results  are  in  agreement  with  previous  calculations  with  the  L2 
norm  error  between  the  channel  packet  method  without  absorbing  boundary  condi¬ 
tions  and  the  channel  packet  method  with  absorbing  boundary  conditions  is  on  the 
order  of  6  x  10“^.  This  is  a  small  error  in  comparison  to  the  dramatic  reduction  in 
computational  effort  the  combination  of  the  channel  packet  method  with  absorbing 
boundary  conditions  produces.  The  grid  for  computing  the  correlation  function  was 
reduced  by  a  factor  of  16.  from  256  x  256  to  64  x  64.  Overall,  the  result  of  the  grid 
reduction  was  a  significant  order  of  magnitude  improvement  in  the  time  necessary 
to  compute  the  correlation  function."*^ 

The  collinear  H  +  H2  reaction  was  also  used  to  demonstrate  the  convergence 
of  the  channel  packet  method  with  absorbing  boundary  conditions.  Figure  4.11 
illustrates  the  convergence  of  channel  packet  method  with  absorbing  boundary  con¬ 
ditions’  solution  to  the  accepted  solution  based  on  previous  work.^^’^^'^"’’^^  Again, 
the  order  of  convergence  for  the  channel  packet  method  with  absorbing  boundary 
conditions  is  tested  for  changes  in  the  temporal  resolution.  As  in  the  one  dimensional 
case,  the  order  of  convergence  is  4.  This  is  illustrated  in  Figure  4.12. 

A  simple  system  of  two  coupled  Morse  oscillators  is  used  to  investigate  the 
effects  of  potential  well  depth  and  reaction  masses  on  quantum  reactive  scattering 
for  two  dimensional  potentials  where  there  is  a  barrierless  well  in  the  interaction 
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region.  However,  S-niatrix  elements  for  the  light-light-light,  medium-light-meclium 
and  heavy-light-heavy  mass  configurations  suffer  from  non-  trivial  errors  introduced 
by  absorbing  boundary  condition  reflection.  While  qualitative  information  about 
the  S-matrix  elements  may  be  gleaned,  the  exact  nature  of  the  S-matrix  elements 
remains  obscured  by  the  error  introduced  by  absorbing  boundary  conditions  reflec¬ 
tion. 

Unlike  the  previous  three  mass  configurations,  the  light-heavy-light  mass  con¬ 
figuration  of  the  model  two  coupled  Morse  oscillator  potential  does  not  suffer  from 
significant  absorbing  boundary  conditions  reflection  error.  The  exact  relationship 
between  the  kinetic,  potential,  and  absorbing  boundary  conditions  terms  in  the 
Hamiltonian  that  leads  to  only  negligible  absorbing  boundary  condition  reflection 
is  not  clear.  However,  the  influence  of  well  depth  and  kinetic  coupling  on  quantum 
reactive  scattering  can  be  investigated  using  the  light-heavy-light  configuration  since 
it  does  not  suffer  from  absorbing  boundary  condition  reflection  error.  Increasing  well 
depth  leads  to  a  greater  number  of  and  more  complex  resonances  in  the  S-matrix 
elements.  However,  there  is  no  clear  analytic  connection  between  the  virtual  states 
of  the  harmonic  limit  of  the  two  coupled  Morse  oscillators  and  the  S-matrix  element 
resonances. 

7.1  Further  research 

There  is  a  great  deal  of  research  on  two  dimensional  quantum  reactive 
scattering  that  remains  to  be  done.  The  work  presented  here  provides  a  firm  foun¬ 
dation  for  further  research.  As  we  saw  in  the  6ne  dimensional  case,  potential  wells 
with  barriers  tend  to  trap  long-lived,  quasi-bound  states.  These  states  may  take 
quite  some  time  to  completely  exit  the  interaction  region.  Methods  that  do  not  take 
advantage  of  absorbing  boundary  conditions  would  require  excessively  large  grids 
to  support  the  evolving  wave  packet  until  the  quasi-bound  state  completely  exits 
the  interaction  region.  In  one  dimension,  the  application  of  absorbing  boundary 
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conditions  eliminated  the  need  for  such  a  large  grid.  However,  this  has  yet  to  be 
demonstrated  for  two  dimensions. 

There  is  also  the  research  into  which  potentials  the  channel  packet  method  with 
absorbing  boundary  conditions  is  applicable.  The  combined  method  apparently  fails 
to  work  for  potential  surfaces  where  there  is  a  well  in  the  interaction  region  and  the 
scattering  atoms  are  in  anything  other  than  a  light-heavy-light  mass  configuration. 
It  appears  that  reflection  from  the  absorbing  boundary  conditions  are  the  culprit 
though  the  mixing  of  eigenstates  of  the  full  Hamiltonian  as  in  Equation  6.10. 

Finally,  there  is  the  extension  of  the  channel  packet  method  with  absorbing 
boundary  conditions  to  scatter  in  three  dimensions.  Current  research  in  three  di¬ 
mensional  reactive  scattering  is  limited  to  very  simple  reactions  like  H  +  H2.  While 
computational  power  continues  to  increase,  the  current  methods  used  in  three  di¬ 
mensional  scattering  again  rely  on  large  grids  which  are  inefficient. 
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Appendix  A.  Derivation  of  the  formula  for  scattering  matrix 

elements 

The  formula  in  the  channel  packet  method  for  scattering  matrix  elements 
given  in  Equation  2.30  is  derived  from  the  orthogonality  relationship 


=  {k'y,i'\s\k,n) 


(A.l) 


Inverting  the  final  line  in  Equation  A. 2  shows  that  all  that  is  necessary  to  compute 
S-matrix  elements  is  to  evaluate  the  scalar  product 


7+)-  (A.2) 

Evaluating  Equation  A. 2  begins  by  defining  the  state  |A£+)  as  the  Fourier 
transform  of  the  time  evolution  of  the  reactant  Mpller  state.  The  state  |A^+)  is 
given  by 

\^E+)  =  j  dfexp  (^i^^exp  (A.3) 

Substituting  the  expansion  of  in  terms  of  |A:-y,7)  given  in  Equation  2.22  and 
expanding  H  in  terms  of  T  and  V  gives 

\^E+)  =  {k^)  dt  exp  (^-i  ^  )  ^xp  ■  (A.4) 

Now  consider  the  second  integral  in  Equation  A.4.  The  integral  of  two  complex 
exponentials  yields  a  delta  function  given  by 


A-1 


/"”rf(exp(-i 


exp  I  i 


where  the  substitution 


has  been  used. 


T-y  =  -TT^ 

2fi~y 


(A.6) 


Two  more  properties  of  delta  functions  are  used  to  simplify  Equation  A. 5. 


S  1 -n  1  =  (n) 


(A.7) 


=5^|.5(a-6)  +  <S(a  +  6)l 


(A.8) 


reduce  the  complexities  in  Equation  A. 5  to 


■Ty  +  Vy\  (  ■Et\  'iTTfJ..^ 


[<J  +  5  (^'9) 


The  two  delta  functions  in  Equation  A. 9  reduce  the  integral  over  dk.y  in  Equation 
A. 4  to  a  simple  sum  given  by 


l^£;+) 


2TTfi-y 


[T}+{+ky)  \+k^n+)  +  rj+{-k^)  |-/:^,7+)], 


(A.IO) 


where  the  ±  in  front  of  the  k-y  explicitly  labels  positive  (+)  and  negative  (— )  mo¬ 


menta. 


Evaluating  the  scalar  product  where  Equation  2.22  is  used  to  ex¬ 

pand  in  terms  of  |fcy,7')  to  obtain  S-matrix  elements  leads  to 

(^y)  7+ (+^7)  (^y’V-|+^%>7+) 

+??!.  (/jy)  'n+  (“^7)  (^y57^“l“^757+)]-  (A. 11) 

The  orthogonality  relationship  in  Equation  2.24  reduces  the  integral  in  Equation  A. 11 
to 


{r_'|/ix(£)> 


(+^y)  (“^7)  ^-k-, 


(~^y)  'n+  ("■^7) 


f 


The  remaining  derivation  follows  Equation  2.28  on  page  2-10. 


(A.12) 
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Appendix  B.  Jacobi-to-bond  transformation  for  momentum 

In  the  transformation  from  the  Jacobi  representation  to  the  bond  repre¬ 
sentation,  the  momenta  are  no  longer  separable.  The  transformation  from  Jacobi 
coordinates  to  bond  coordinates  is  given  by 


(B.l) 


where  <j)  is  the  mass  related  factor  given  by 


4>  = 


rric 
nib  + 


(B.2) 


This  transformation  is  for  the  A  +  BC  channel  where  is  the  distance  between  the 
masses  in  the  diatom  BC  and  Ra  is  the  distance  between  the  atom  A  and  the  center- 
of-mass  of  the  diatom.  The  derivation  of  the  momentum  transformation  is  also  for 
the  A  -b  BC  channel.  For  the  AB  +  C  channel,  the  derivations  and  transformations 
are  similar. 

For  calculating  the  momentum  transformation,  the  starting  point  is  the  La- 
grangian  formulation 

L  =  T-  V;  (B.3) 

where  L  is  the  Lagrangian,  T  is  the  kinetic  energy  and  V  is  the  potential  energy. 
For  the  quantum  reactive  scattering  studied  in  this  research  project,  the  potential 
V  is  independent  of  the  velocity  and  it  is  omitted  from  the  rest  of  the  discussion  of 
the  transformation.  The  conjugate  momenta  P,  is  derived  from  the  Lagrangian  via 


(B.4) 
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where  cj,'  is  the  time  derivativ'e  of  the  coordinate  cj.  In  Jacobi  coordinates,  the 
Lagrangian  for  a  free  particle 


'2flr 


where  is  the  mass  of  atom  a,  /Xr  is  the  reduced  mass  given  by 


m^mc 

l^r  —  j  ? 
rrth  +  rUc 

and  Pr  and  Pr  are  the  non-relativistic  momenta  given  by 


Pr  =  rUaR  and  Pr  =  Hri'- 


(B.5) 


(B.6) 


(B.7) 


Combining  together  equations  B.5  and  B.7  gives 

L  =  (B.8) 

Now.  using  the  coordinate  transformations  from  equation  B.2  together  with  equa¬ 
tion  B.8,  yields 

L  =  im,  {X  P  ci>Yf  +  ^-^irY\  (B.9) 

Using  the  relationship  in  equation  B.4,  the  conjugate  momenta  Px  and  Py  are  given 

by 

Px  =  |^  =  im.(2A:+2ey)  (B.IO) 

Py  =  |^=lm.(2^X+2y)  +  ^.F.  (B.ll) 
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Combining  terms  together  with  the  transformation  in  equation  B.2  yields  the  final 
transformation, 


Px  —  rriaR  (B.12) 

Py  =  ma<t)R  +  r  [ilia  +  1-i.r  —  ■  (B.13) 

In  the  A-\-  BC  channel,  Px  may  be  interpreted  as  the  relative  momentum.  However, 
the  interpretation  of  Py  is  clouded  by  mixing  between  relative  and  internal  momenta. 
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Appendix  C.  Code  overview 

The  code  implementing  the  application  of  the  channel  packet  method  to¬ 
gether  with  absorbing  boundary  conditions  to  quantum  scattering  in  one  and  two 
dimensions  is  briefly  overviewed  here.  The  code  is  written  in  FORTRAN  ( i  with 
extensions  made  available  by  the  Silicon  Graphics  compiler. 

The  fast  Fourier  transform  (FFT)  used  throughout  this  code  is  the  Temperton 
FFT  written  by  the  CECAM  group.^®  The  two  dimensional  front  end  to  the  FFT 
was  written  by  the  author.  The  original  two-dimensional  code  was  written  by  Dr. 
David  Weeks  during  a  post-doctoral  fellowship  at  Notre  Dame  University  under  the 
direction  of  Dr.  David  Tannor  and  was  altered  to  include  the  application  of  absorbing 
boundary  conditions  by  the  author.  Both  the  one  and  two  dimensional  codes  rely 
upon  the  use  of  individual  modules  of  code  compiled  into  a  single  executable  using 
the  Unix  make  command.  The  initialization  of  the  code  variables  is  done  using  a 
header  file  and  all  variables  are  common.  No  screen  input  is  required  as  the  input 
data  is  contained  in  a  data  file  read  into  the  program  at  the  start. 

C.l  One  dimensional  code 

The  code  implementing  the  application  of  the  channel  packet  method  to¬ 
gether  with  absorbing  boundary  conditions  to  quantum  scattering  in  one  dimension 
was  written  by  the  author.  The  module  titles  are  in  boldface  with  a  brief  explana¬ 
tion  of  its  function. 

iDFFT.f:  This  module  contains  the  Temperton  FFT  and  the  initialization 
modules  for  the  FFT. 

iDScatter.h:  Header  file  containing  all  the  variables  used  in  the  program. 
All  the  variable  are  common. 

iDScatter.f:  Main  program  file  which  calls  each  of  the  modules  in  turn. 
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iDVariables.f:  This  module  reads  in  the  input  data  file  containing  the  grid 
and  wave  packet  initialization  data. 

iDGrid.f:  This  module  creates  the  coordinate  and  momentum  representation 

grids. 

iDWave.f:  This  module  computes  the  initial  wave  packets. 

iDPotential.f:  This  module  creates  the  potential  array  based  on  the  coordi¬ 
nate  grid. 

iDMatrix.f:  This  module  implements  the  split-operator  array  representation 
of  the  full  scattering  Hamiltonian  using  the  coordinate  and  momentum  grids  together 
with  the  potential  array.  It  is  here  that  absorbing  boundary  conditions  are  introduced 
into  the  Hamiltonian. 

iDNuscreen.f  and  iDGraphs.f:  The  first  module  initializes  the  graphics 
screen  and  the  second  module  draws  to  the  screen  using  the  OpenGL  libraries  on 
a  Silicon  Graphics  workstation.  The  modules  are  specific  to  the  Silicon  Graphics 
operating  system  and  are  not  for  use  on  any  other  platform. 

iDReactant.f:  This  module  computes  the  reactant  Mpller  state  in  the  manner 
outlined  in  section  3.1.2.  On  a  Silicon  Graphics  workstation,  the  evolving  wave  packet 
can  be  displayed  on  the  screen  through  the  IDGraphs.f  module. 

IDProduct.f:  This  module  computes  the  product  Mpller  state  in  the  manner 
outlined  in  section  3.1.2.  On  a  Silicon  Graphics  workstation,  the  evolving  wave 
packet  can  be  displayed  on  the  screen  through  the  IDGraphs.f  module. 

iDNuWay.f:  This  module  propagates  the  product  Mpller  state  forw'ards  and 
backwards  in  time  calculating  the  correlation  function  at  every  step.  On  a  Sili¬ 
con  Graphics  workstation,  the  evolving  wave  packet  can  be  displayed  on  the  screen 
through  the  IDGraphs.f  module. 

IDSmatrix.f:  This  module  computes  the  S-matrix  elements  according  to 
equation  2.29.  The  module  iDFTCR.f  actually  computes  the  discrete  Fourier  trans- 


form  of  the  correlation  function.  The  resulting  S-matrix  elements  are  written  to  a 
file  and  the  program  terminates. 

C.2  Two  dimensional  code 

The  calculation  of  S-matrix  elements  for  two  dimensional  reactive  quan¬ 
tum  scattering  is  actually  done  using  three  separate  codes:  Mpller  state  calculation, 
correlation  function  calculation  and  S-matrix  calculation.  The  code  computing  the 
Mpller  states  remains  unaltered  from  Dr.  Weeks’  original  work.  The  code  computing 
the  correlation  function  is  based  on  Dr.  Weeks’  original  work  but  altered  to  include 
absorbing  boundary  conditions.  The  code  computing  S-matrix  elements  was  written 
by  the  author  for  this  project. 

C.2.1  M0ller  state  calculation  code.  The  module  titles  are  in  boldface 

with  a  brief  explanation  of  its  function. 

TemptFFT.f:  This  module  contains  the  Temperton  FFT  and  the  initializa¬ 
tion  modules  for  the  FFT. 

2DScatter.h:  Header  file  containing  all  the  variables  used  in  the  program. 
All  the  variable  are  common. 

2DScatter.f:  Main  program  file  which  calls  each  of  the  modules  in  turn. 

2DInitVar.f:  This  module  reads  in  the  input  data  file  containing  the  grid  and 
wave  packet  initialization  data. 

2DInitGrid.f:  This  module  creates  the  coordinate  and  momentum  represen¬ 
tation  grids. 

2DPotential.f:  This  module  creates  the  potential  array  based  on  the  coordi¬ 
nate  grid.  This  module  calls  another  module  that  contains  the  code  for  the  potential 
energy  surface.  The  collinear  H  +  H2  potential  energy  surface  is  contained  in  the 
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module  hh2pot.f  and  the  two  coupled  Morse  oscillator  potential  energy  surface  is 
contained  in  the  module  morsepot.f. 

2DInitWave.f;  This  module  computes  the  initial  wave  packets  which  is  more 
complicated  in  the  2D  case.  At  each  grid  point,  depending  on  the  channel,  the 
coordinates  are  transformed  into  the  appropriate  set  of  Jacobi  coordinates.  The  Ja¬ 
cobi  coordinates  are  then  passed  to  a  module  which  computes  the  internal  vibrational 
eigenstate.  For  the  collinear  H+  H2  potential  energy  surface,  the  internal  vibrational 
eigenstates  are  not  analytic.  The  internal  Hamiltonian  is  represented  in  a  basis  of 
Morse  oscillator  functions  and  then  diagonalized.  The  internal  vibrational  eigen¬ 
state  is  computed  from  a  linear  combination  of  Morse  oscillator  functions  according 
to  the  transformation  matrix  which  diagonalized  the  internal  Hamiltonian.  For  the 
two  coupled  Morse  oscillator  potential  energy  surface,  the  internal  vibrational  eigen¬ 
states  are  known  analytically  and  the  matrix  representation  and  diagonalization  are 
omitted.  In  either  case,  it  is  important  to  note  that  the  code  numbers  the  internal 
eigenstates  starting  at  1  and  not  0.  The  initial  wave  packet  is  constructed  from 
a  direct  product  of  the  internal  vibrational  eigenstate  and  a  linear  combination  of 
eigenstates  of  the  relative  Hamiltonian  in  Jacobi  coordinates.  It  is  then  propagated 
analytically  to  ±r  and  then  transformed  into  bond  coordinates. 

2DMatrix.f;  This  module  implements  the  split-operator  array  representation 
of  the  full  scattering  Hamiltonian  using  the  coordinate  and  momentum  grids  together 
with  the  potential  energy  surface. 

2DSplitOp.f:  This  module  computes  the  Mpller  state  from  the  initial  wave 
packet  using  the  split  operator  approach.  At  the  end  of  the  propagation,  the  Mpller 
state  is  written  to  a  file  and  the  program  terminates. 

C.2.2  Correlation  function  calculation  code.  The  module  titles  are 

in  boldface  with  a  brief  explanation  of  its  function.  The  overviews  of  the  mod- 
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ules  TemptFFT.f  2DScatter.h,  2DScatter.f,  2DInitVar.f,  2DInitGrid.f  and 
2DPotential.f  remain  unchanged  from  the  previous  section. 

2DInitWave.f:  This  module  simply  reads  in  the  reactant  and  product  Moller 
states. 

2DMatrix.f:  This  module  implements  the  split-operator  array  representation 
of  the  full  scattering  Hamiltonian  using  the  coordinate  and  momentum  grids  together 
with  the  potential  energy  surface.  This  code  has  been  altered  to  include  absorbing 
boundary  conditions. 

2DSplitOp.f:  This  module  propagates  the  reactant  Mpller  state  forwards 
and  backwards  in  time  using  the  split  operator  approach.  At  each  time  step,  the 
correlation  function  is  computed  and  stored.  At  the  end  of  the  propagations,  the 
correlation  function  is  written  to  file  and  the  program  terminates. 

C.2.3  S-matrix  calculation  code.  The  module  titles  are  in  boldface 

with  a  brief  explanation  of  its  function. 

2DScatter.h:  Header  file  containing  all  the  variables  used  in  the  program. 
All  the  variable  are  common. 

2DSmatrix.f:  This  is  the  main  module.  It  first  calls  2DInitVar.f  to  input 
all  the  necessary  information.  The  main  body  of  code  computes  the  expansion 
coefficients  necessary  to  normalize  the  S-matrix  elements,  computes  the  discrete 
Fourier  transform  of  the  correlation  function  and  calculates  and  writes  to  a  file  the 
S-matrix  elements. 

2DInitVar.f:  This  module  not  only  reads  in  the  initial  grid  and  wave  packet 
conditions,  it  also  reads  in  the  correlation  function. 
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Appendix  D  S-matrix  elements  for  two  coupled  Morse  oscillators  in  a 


light-heavy-light  mass  configuration 


The  following  16  plots  illustrate  the  probability  of  reaction  for  two  coupled  Morse 
oscillators  in  a  light-heavy-light  mass  configuration  where  the  well  depth,  i.e. 
dissociation  energy,  is  changed  from  0.1 1  au  to  0.26  au  in  steps  of  0.01  au. 


Figure  D.l  Absolute  value  of  the  correlation  function  for  the  two  coupled  Morse 
oscillator  LHL  mass  configuration  for  the  reaction 
A  +  BC(v  =  0) -»AB(v’=0)-i-C  where  the  disassociation  energy  is 

0.11  au. 


Figure  D.2  Absolute  value  of  the  correlation  function  for  the  two  coupled  Morse 
oscillator  LHL  mass  configuration  for  the  reaction 
A+BC(v  =  0) AB(v’=0)+C  where  the  disassociation  energy  is 

0.12  au. 
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Relative  Translational  Energy  (eV) 


D.3  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  LHL  mass 
configuration  for  the  reaction  A  +  BC{v  =  0)  ^  AB(v’=  0)  +  C  where  the 

disassociation  energy  is  0.13  au. 


Figure  D.4  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  LHL  mass 
configuration  for  the  reaction  A  +BC(v  =  O)  — ^  AB{y'=  O) +C  where  the 

disassociation  energy  is  0.14  au. 
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Figure  D.5  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  LHL  mass 
configuration  for  the  reaction  A  +  BC{y  =  O)  — >  AB(v’=  0) +C  where  the 

disassociation  energy  is  0.15  au. 


Figure  D.6  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  LHL  mass 
configuration  for  the  reaction  A  +  BC{v  =  0)  ->  AB{v’=  0)  +C  where  the 

disassociation  energy  is  0. 16  au. 
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disassociation  energy  is  0.17  au. 


Figure  D.8  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  LHL  mass 
configuration  for  the  reaction  A  +  BC{y  =  0)  — >  0)  +C  where  the 

disassociation  energy  is  0.18  au. 
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Figure  D.9  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  LHL  mass 
configuration  for  the  reaction  A  +BC(v  =  0)  — >  AB[v*=  0)  +C  where  the 

disassociation  energy  is  0. 19  au. 


Figure  D.  10  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  LHL  mass 
configuration  for  the  reaction  A  +BC(v  =  O)  — >  AB(v’=  0) +C  where  the 

disassociation  energy  is  0.20  au. 
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FigureD.il  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  LHL  mass 
configuration  for  the  reaction  A  +  BC{y  =  O)  ^  AB{y=  O) +C  where  the 

disassociation  energy  is  0.21  au. 


Relative  Translational  Energy  (eV) 

Figure  D.  12  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  LHL  mass 
configuration  for  the  reaction  A  +  BC(v  =  O)  — >  AB(v'’=  0)  +C  where  the 

disassociation  energy  is  0.22  au. 
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Figure  D.  13  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  LHL  mass 
configuration  for  the  reaction  A  +  BC[v  =  0)  — >  AB(^v’=  0)  +C  where  the 

disassociation  energy  is  0.23  au. 


Figure  D.  14  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  LHL  mass 
configuration  for  the  reaction  A  +  BC{y  =  0)  — >  AB{y'—  0)  +C  where  the 

disassociation  energy  is  0.24  au. 
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Figure  D.  15  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  LHL  mass 

configuration  for  the  reaction  A  +  BC(w  =  0)  AB(v’=  0)  +C  where  the 

disassociation  energy  is  0.25  au. 


Figure  D.  16  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  LHL  mass 
configuration  for  the  reaction  A  +  BC{y  =  0)  — >  AB{y*=  0)  +C  where  the 

disassociation  energy  is  0.26  au. 
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The  following  11  plots  illustrate  the  probability  of  reaction  for  two  coupled 
Morse  oscillators  in  a  light-heavy- light  mass  configuration  where  the  kinetic  energy 
coupling  constant  is  changed  from— 1  (fully  coupled)  to  0  (uncoupled)  in  steps  of  0.1. 


Figure  D.17  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  LHL  mass 
configuration  for  the  reaction  A  +  BC{v  =  0)  — >  AB{v^=  0)  H-C  where  the 
disassociation  energy  is  0.20  au  and  the  kinetic  energy  coupling  constant 
is  -1. 


Relative  Translational  Energy  (eV) 


Figure  D.  18  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  LHL  mass 
configuration  for  the  reaction  A  +BC{v  =  0)  -4  AB(w’=  0) +C  where  the 
disassociation  energy  is  0.20  au  and  the  kinetic  energy  coupling  constant 
is  -0.9. 
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Figure  D.  19  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  LHL  mass 
configuration  for  the  reaction  A+BC{v  =  0)  ^  AB{v'=  0)  +C  where  the 

disassociation  energy  is  0.20  au  and  the  kinetic  energy  coupling  constant  is 

-0.8. 


Figure  D.20  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  LHL  mass 
configuration  for  the  reaction  A  +  BC{v  =  0)  — >  AB{v’=  0) +C  where  the 

disassociation  energy  is  0.20  au  and  the  kinetic  energy  coupling  constant  is 
-0.7 
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D.21  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  LHL  mass 
configuration  for  the  reaction  A  +  jBC(v  =  O)  — >  AJ5(v’=  0) +C  where  the 
disassociation  energy  is  0.20  au  and  the  kinetic  energy  coupling  constant  is 
-0.6. 
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Figure  D.22  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  LHL  mass 
configuration  for  the  reaction  A  +  BC(v  =  0)  — >  AJB(y  =  0)  +C  where  the 
disassociation  energy  is  0.20  au  and  the  kinetic  energy  coupling  constant  is 
-0.5. 
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Figure  D.24  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  LHL  mass 
configuration  for  the  reaction  A  =  0)  AB(v^=  0) +C  where  the 

disassociation  energy  is  0.20  au  and  the  kinetic  energy  coupling  constant  is 
-0.3. 
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D.25  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  LHL  mass 
configuration  for  the  reaction  A  +  BC{v  =  0)  -4  AB(v’=  0) +C  where  the 
disassociation  energy  is  0.20  au  and  the  kinetic  energy  coupling  constant  is 
-0.2. 


Figure  D.26  Probability  of  reaction  for  the  two  coupled  Morse  oscillator  LHL  mass 
configuration  for  the  reaction  A  +  BC(v  =  0)  — >  A.B(v  =  0)  +C  where  the 
disassociation  energy  is  0.20  au  and  the  kinetic  energy  coupling  constant  is 
-0.1. 
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